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ABSTRACT 

Engineers are challenged to produce better designs in less time and for less cost. 
Hence, to investigate novel and revolutionary design concepts, accurate, high-fidelity 
data must be assimilated rapidly into the design, analysis and simulation process. This 
data assimilation should consider diverse mathematical modeling and multi-discipline 
interactions necessitated by concepts exploiting advanced materials and structures. 
Integrated high-fidelity methods with diverse engineering applications provide the 
enabling technologies to assimilate these high-fidelity, multi-disciplinary data rapidly at 
an early stage in the design. These integrated methods must be multifunctional, 
collaborative and applicable to the general field of engineering science and mechanics. 

Multifunctional methodologies and analysis procedures are formulated for 
interfacing diverse domain idealizations including multi-fidelity modeling methods and 
multi-discipline analysis methods. These methods, based on the method of weighted 
residuals, ensure accurate compatibility of primary and secondary variables across the 
domain interfaces. Methods are developed for scalar-field and vector-field problems in 
engineering science with extensions to multidisciplinary problems. Results are presented 
for the scalar- and vector-field developments using example patch test problems. In 
addition, results for torsion, thermal, and potential flow problems are presented to 
demonstrate further the effectiveness of the scalar-field development. Results for plane 
stress and plane flow problems are presented for the vector-field development. Results 
for all problems presented are in overall good agreement with the exact analytical 


solution or the reference numerical solution. 



IV 


The multifunctional methodology presented provides an effective mechanism by 
which domains with diverse idealizations are interfaced. This capability rapidly provides 
the high-fidelity data needed in the early design phase. Moreover, the capability is 
applicable to the general field of engineering science and mechanics. Hence, it provides 
a collaborative capability that accounts for interactions among engineering analysis 
methods. 
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CHAPTER I 
INTRODUCTION 

1.1. MOTIVATION 

The analysis of revolutionary aerospace and ground vehicles relies heavily on 
accurate, efficient and robust computational methodologies such as the finite element and 
finite difference methods. To investigate novel and revolutionary design concepts, 
accurate, high-fidelity data must be assimilated rapidly into the design, analysis and 
simulation process. This data assimilation should consider mathematical modeling 
approximations ranging from simple handbook equations, empirically derived relations, 
spreadsheets, and design charts to complex continuous and discrete simulation models. 

In addition, the data assimilation needs to consider associated multi-discipline 
interactions necessitated by advanced design concepts exploiting multifunctional 
materials and leading to multifunctional structures. Rapid discipline-centric modeling 
techniques allow high-fidelity design trades between cost and performance, and based on 
the insight provided by these simulations, design uncertainties and risk assessment may 
be evaluated. Integrated multi-discipline analyses allow the assessment of the effects of 
multidisciplinary coupling on the system response. New computing systems and 
alternative computing strategies have presented new opportunities for optimal design, 
analysis, and simulation of aerospace systems. However, integrated high-fidelity 
methods with diverse engineering applications provide the enabling technologies to 
assimilate high-fidelity, multi-disciplinary data rapidly at an early stage in the design. 


The journal model for this dissertation is the AIAA Journal. 
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These integrated methods must be multifunctional, collaborative and applicable to the 
general field of engineering science and mechanics. 

To understand the impact of these integrated methods, the three concomitant 
attributes, namely, multifunctional, collaborative, and engineering science and 
mechanics, must be described. In the context of this work, multifunctional 
characterization has been adopted from the description of new and innovative materials 
and structures with multiple capabilities. These systems, referred to as multifunctional 
materials and structures, respectively, have several desirable simultaneous properties and 
many diverse disciplinary applications. The systems will adapt, react and evolve in 
changing environments, and their use will result in a combined system with enhanced 
capabilities at less cost and weight. Likewise, multifunctional methods refer to 
computational methodologies that have multiple capabilities such as multiple fidelity 
modeling, multiple approximation analysis and multidisciplinary analysis. The methods 
are computationally efficient while preserving solution accuracy and are applicable to a 
wide range of applications in engineering science. Their use in the combined analysis of 
complex configurations promises to provide enhanced computational and engineering 
capability at less cost and in less time. With these attributes, a multifunctional method 
may address the diverse modeling and analysis needs of evolving systems perhaps using 
a hierarchical approach including error analysis and risk assessment. 

The collaborative aspect of the computational methods provides a mechanism by 
which two or more physical domains are integrated or interfaced and by which two or 
more methods or algorithms are shared or interfaced. It is through this interfacing that 
the diverse attributes create a unified framework that far exceeds the capability of an 
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individual method. Collaborative methods may integrate domains of different 
discretization fidelity, analysis approximations, or disciplines. An example of a 
collaborative method is adaptive dynamic relaxation. Explicit direct time integration 
algorithms are well-known for their computational efficient, low-memory requirements, 
low computational cost per solution step and direct mapping to massively parallel 
processing (MPP) systems. Adaptive dynamic relaxation techniques exploit these 
features to determine the quasi-static or steady-state response of a structure without 
relying on traditional methods requiring the solution of the large sparse matrices. 
Collaborative methods provide a mechanism by which the aggregate cost savings related 
to computational and modeling requirements are reduced, and analyses, previously 
intractable, may be performed. As in the case of the multifunctional materials or 
structures, these methods adapt, react and evolve in the changing environments of 
engineering science. Engineering science covers the broad perspective of engineering 
and includes the integrated application of engineering principles, science, mathematics, 
numerical analysis and non-deterministic methods. Problems in fluid flow, solid 
mechanics, thermal analysis, and constitutive modeling are representative of those in 
engineering science. Engineering science has a multidisciplinary emphasis, and future 
methods applicable to the field should possess multifunctional characteristics and a 
collaborative nature to further enhance their analysis capabilities and to advance the state- 
of-the-art in engineering design. 

Multifunctional collaborative methods should address four typical steps of 
analysis and design, namely, (1) representation or modeling of the geometry, (2) 
knowledge-based selection and development of appropriate mathematical models ( i.e ., 
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idealization/discretization), (3) solution of the mathematical model (continuous and/or 
discrete), and (4) interrogation/assessment of the results. These steps are briefly outlined 
in Appendix A. Methodology and analysis procedures that address these basic steps 
provide the foundation for enhanced integrated design and analysis tools within the realm 
of engineering science. Such multifunctional methodology should allow interaction 
between and collaboration with the analyst and designer, among different mathematical 
modeling approximations of the physical phenomena, and among multiple engineering 
disciplines. A major feature of the methodology is the transfer of data across the 
respective interface, whether the interface is one among diverse mathematical 
approximations or among diverse disciplines. Computational issues associated with 
individual modeling approaches and disciplines are magnified in number and significance 
due to the intricate couplings manifesting themselves as a by-product of their interfacing. 

Multi-fidelity modeling approaches provide benefits in all of the major steps of 
analysis and simulation. These approaches are often characterized by the use of different 
approximations among multiple domains of the same continua and multiple domains 
involving different continua ( e.g . , fluid-structure interaction). Analytical and closed- 
form solutions for specific geometries and configurations are often used to eliminate 
constraints placed on the analysis due to geometry considerations. Rapid modeling 
approaches facilitate the discretization of geometry by providing a capability to model 
regions of interest, independently, increasing the discretization fidelity or enhancing the 
mathematical approximation only in the desired domains. Thus, for multi-fidelity finite 
element modeling approaches, complex and often unsuitable mesh transitioning, 
generated manually or using automatic mesh generators, is limited. In addition, multi- 
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fidelity approaches have been developed that allow for the discretization of parts or 
components across geographically dispersed locations with minimal concern for the 
discretization of the parts along common boundaries or interfaces. Additional research 
has provided for accommodation of slight anomalies in the geometric representation 
provided by the independently discretized parts as well as parametric definition of the 
interface geometry between parts. Multi-fidelity modeling approaches benefit the 
solution of the discretized system in that the system size using a multi-domain approach 
for global/local modeling may be smaller for a given level of solution accuracy than the 
system obtained by standard practices. In addition, in component modeling, the 
associated matrices may be reduced by static condensation, which reduces the size and 
subsequent solution time of the overall system of equations. Multi-fidelity modeling 
approaches allow for the visualization and interrogation of the results only in regions of 
interest. Post-processing of secondary results such as stresses and failure parameters may 
be isolated to these regions and dynamically computed as the need arises. By reducing 
the modeling, computational and visualization time of simulations of aerospace 
structures, multi-fidelity modeling approaches promise to enhance the viability of high- 
fidelity analyses early in the design process. 

Multidisciplinary coupling approaches involve the interfacing of different 
disciplines to account for their interactions and impact on the overall system response. 
There are myriad approaches, for example, any combination of approaches that couple 
the fluids, thermal, structures, and acoustic disciplines. The traditional independent 
approach for multidisciplinary analysis involves loosely coupling the disciplines through 
sequential execution of single discipline analyses. Typically this approach requires 
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several iterations among the different analysis methods and analysts and is relatively 
inefficient because the discipline specific models are generally incompatible and require 
extensive post-processing after each single discipline analysis to transfer (or interface) 
data to the next analysis model. Aeroelastic analysis as an interdisciplinary problem, 
requires the coupling of the aerodynamic and structural responses. The use of different 
spatial discretization procedures and potentially different mathematical modeling 
approximations for the aerodynamic model and the structures model gives rise to the 
interfacing problem of transferring computed data between the two grid systems. 
Moreover, the same issues are prevalent in fluid-thermal-structural analyses and 
structural-acoustic analyses. Suitable methodology for addressing these types of 
interfacing problems has been developed by many researchers. 

The overarching purpose of this research is to investigate multifunctional 
collaborative methods, as described herein, that address the engineering design and 
analysis needs of multidisciplinary problems in engineering science. This research 
focuses on the fundamental relationships among underlying engineering science and 
mechanics principles, computational methods and multi-fidelity models, and methods 
using basic problems from continuum mechanics. Given its broad applicability with 
respect to the field of engineering science, continuum mechanics forms the foundation for 
the multifunctional collaborative methods developed in this work. Hence, for 
completeness and to establish notation, basic concepts of continuum mechanics are 
presented briefly in the next section. 
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1.2. CONTINUUM MECHANICS FOUNDATIONS 

Continuum mechanics is the branch of physical sciences concerned with the 
deformations and motions of continuous material media under the influence of external 
effects 1 . The effects that influence the bodies appear in the form of forces, 
displacements, and velocities which arise from contact with other bodies, gravitational 
forces, thermal changes, chemical interactions, electromagnetic effects, and other 
environmental changes. In this work, bodies subject to forces of mechanical origin 
and/or thermal changes are of primary concern. General principles in the form of integral 
or differential equations govern the deformation and motion of the continuum. Hence, 
approximation methods and associated concepts are introduced in addition to the basic 
concepts of continuum mechanics. 

1.2.1. General Principles of Continuous Media 

A medium can be generally categorized as a fluid or a solid. A fluid can be 
loosely defined as a continuum that does not require external forces to maintain its 
deformed shape. When highly compressible it is called a gas and when essentially 
incompressible, it is called a liquid. A solid can be loosely defined as a continuum that 
requires external forces to maintain its deformed shape. According to its behavior, a 
solid may be called elastic, plastic, viscoelastic, thermoelastic, etc. Usually it is assumed 
to have a uniform density . When a medium deforms, the small volumetric elements 
change position by moving along space curves. Their positions as functions of time can 
be specified either by the Lagrangian ( X t = Xf(x(,t) for i= 1, 2, 3) or Eulerian 

description ( x,- = x ; - {X(,t ) ). In the Lagrangian description, each particle is tracked in 


terms of its initial position with respect to a fixed reference system, X t , and time. In the 
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Eulerian description, the motion is expressed in terms of the instantaneous position vector 
with respect to a moving reference system, x t , and time. 

Classical continuum mechanics rests upon equations expressing the balances of 
mass, linear momentum, angular momentum, energy, and entropy in a moving body . 
These balance laws apply to all material bodies, whether fluid or solid in composition, 
and each gives rise to a field equation. These balance laws are as follows: 

i. Principle of conservation of mass 

ii. Principle of conservation of linear momentum 

iii. Principle of conservation of angular momentum 

iv. Principle of conservation of energy 

v. Principle of entropy 

The principle of conservation of mass states that when the total mass of the body is 
unchanged for an arbitrarily small neighborhood of each material point, the mass is 
considered to be conserved locally. The conservation of linear momentum represents 
Newton’s second law and governs the motion of the continuum under the influence of the 
external effects. The principle of conservation of angular momentum is used to show 
symmetry of the stress tensor for many engineering materials, and the stress tensor 
describes the state of stress of the continuum. The principle of conservation of energy, 
also called the first law of thermodynamics, states that energy is conserved if the time 
rate of change of the kinetic and internal energy is equal to the sum of the rate of work of 
the external forces and all the other energies entering or leaving the body. The second 
law of thermodynamics is automatically satisfied and includes the change in entropy of 
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the continuum. More detailed descriptions of these balance laws are presented in Chapter 
III. 


In deriving the governing equations, the starting point is a statement of the 
conservation principle applied to a “control volume” to develop the integral form of the 
equation and extract the differential form by using the divergence theorem. A control 
volume has a fixed volume in space; its boundary does not deform but allows mass 
transfer through it. In contrast, a material volume contains the same quantity of material 
at all times; its boundary can deform, and it does not allow mass transfer. 

As the continuum moves, in general, properties change with time and space. The 
material derivative (substantial or total) must account for these changes depending on the 
method of description used. Consider the scalar property as tp, for the Lagrangian 
description, the material derivative is: 


dXj _ d0 

dt dt dX j dt dt 

For the Eulerian description, the material derivative is: 

d<p(xj,t) _ dtp dtp dx\ dtp dx 2 dtp dx 3 

dt dt dx\ dt d%2 dt dx 3 dt 

dtp dtp 
= — + v,— — 
dt 1 dx{ 


= — + v- V(|) 
dt 


The general conservation equation may be written in integral form or differential 
form in conservative or divergence form. However, when considering the differential 
form, an equivalent representation is often obtained by working out the divergence 
operator and introducing the material derivative. This leads to a non-conservative form 
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of the differential equation. Although the conservative and non-conservative forms of the 
differential equations of the conservation principles are equivalent from a mathematical 
point of view, they will not necessarily remain so when a numerical discretization is 
performed. The general form of the conservation law is said to be written in conservative 
or divergence form. The importance of the conservative form in a numerical scheme lies 
in the fact that, if not properly taken into account, a discretization of the differential 
equations will lead to a numerical scheme in which all the mass fluxes through the mesh - 
cell boundaries will not cancel; hence, the numerical scheme will not keep the total mass 
constant 4 . 

1.2.2. Mathematical Approximations 

Mathematical problems frequently encountered in engineering science may be 
classified as boundary- value and initial -value problems based upon the existence of one 
or more supplementary conditions. The differential equation describes a boundary- value 
problem if the dependent variable and possibly its derivative are required to have 
specified values on the domain boundary. The differential equation describes an initial- 
value problem if the dependent variable and possibly its derivative are specified initially 
(, i.e ., t=0). Initial-value problems are generally time dependent. 

Partial differential equations governing the motion of general continua are often 
of the canonical form Au xx + Bu + Cuyy = 0 where the coefficients A, B, and C are 

real constants, u represents a field variable, and the subscripts, x andy, denote partial 
differentiation with respect to the independent variables, x and y. The character of this 
quasi-linear, second-order, partial differential equation is determined by the sign of the 

discriminant, B -4 AC. The partial differential equation is 
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elliptic 

for B 2 -4 AC < 0 

hyperbolic 

for B 2 -4AC> 0 

parabolic 

for B 2 -4AC = 0 


The full significance of the classification of quasi-linear, second-order partial differential 
equations as elliptic, hyperbolic, or parabolic is beyond the scope of this work. However, 
this classification has proved important for an understanding of the kinds of initial and 
boundary conditions one must furnish along with the partial differential equation in order 
to determine a unique solution. Moreover, solution methods differ markedly from one 
classification to another, which is of particular importance in the field of fluid 
mechanics 6 . For example, boundary conditions are generally imposed all the way around 
a rectangular domain (the x-y region) of a two-dimensional flow when the equation is 
elliptic, and the solution must have no discontinuities in the second derivatives, except 
possibly at singular points where the differential equation is not applicable. Hyperbolic 
and parabolic equations, by contrast, have at least one open boundary; thus, boundary 
conditions are not usually imposed all around the domain under consideration. The 
boundary conditions for at least one variable, usually time, are specified at one end, and 
the system is integrated indefinitely. Certain kinds of discontinuities in the second 
derivatives are admissible across certain curves in such a way that the differential 
equation continues to be applicable in those regions. 

Approximate solutions of differential equations ( e.g ., Ritz, Galerkin, least- 
squares, collocation or in general weighted-residual methods) satisfy only part of the 
conditions of the problem. For example, either the governing equation or the boundary 
conditions may be satisfied only at a few positions rather than at each point. The 
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approximate solution is expanded in a set of known functions with arbitrary parameters. 
Two ways to determine the parameters are the method of weighted residuals (MWR) and 
the variational method. While the MWR and variational methods are only briefly 
discussed here, a more complete discussion of the approaches is given in the literature by 

n 

Finlayson . In MWR, one works directly with the differential equation and boundary 
conditions, whereas in the variational method one tries to satisfy the governing 
differential equation in an average sense using a functional related to the differential 
equations and/or the boundary conditions. MWR encompasses several methods 
(collocation, Galerkin, integral, etc.) and provides a framework to compare and contrast 
methods. Variational methods are not applicable to all problems, and thus suffer a lack 
of generality. MWR is easy to apply whereas variational methods require manipulation 
that can be more complex. 

Variational methods provide a means for the determination of the governing 
equations. In solid mechanics, the principles of virtual work and stationary potential 
energy can be used to derive the governing equations and boundary conditions. The 
principle of virtual work demands that for the state of equilibrium, the work of the 
impressed forces is zero for any infinitesimal variation of the configuration of the system 
that is in harmony with the kinematic constraints. Hence, the variational statement 
implicitly imposes the natural boundary conditions. All work statements are derived 
from classical laws pertaining to the equilibrium of the particle. Moreover, the virtual 
work statement is simply the weak form of the equilibrium equations. For monogenic 
forces, this statement leads to the condition that for equilibrium, the potential energy shall 
be stationary with respect to all kinematically permissible variations. 
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The original differential equation is said to be the strong form of the problem 
while the integral form is typically referred to as the weak form. However, in the strict 
sense, particularly for approximation methods such as the Galerkin method, the weak 
form is obtained by transferring the differentiation from the dependent variable to the test 
functions, which includes the identification of the type of boundary conditions that the 
weak form can admit. The purpose of the transfer of differentiation is to equalize the 
continuity requirements on the dependent variable and the test function. This results in a 
weaker continuity requirement on the solution in the weak form than in the original 
equation. In the process of transferring the differentiation, boundary terms that determine 
the nature of the natural or essential boundary conditions in the solution are obtained. 

The classification of boundary conditions as natural and essential boundary 
conditions plays a crucial role in the derivation of the approximate functions. From 
variational calculus, consider a partial differential equation in the form, 

dF d ( dF "| a 

dll dx dll r dy 
x J x 

where F = F(x,y,u,u x u y ) , u x = dit / dx and u y = du / dy . Transferring the 

differentiation from the dependent variable, u, to the test function, v, yields the weak 
form of the differential equations in the form 

J dF dv dF dv dF ~L , A dF dF V A 
I v ” h ~ ~ — +— — — dxdy-4v- — n x +- — n x ds = 0 
q du dx du x dy du y r du x du y 

It is at this point that the natural and essential boundary conditions are readily identified. 
Generally, specifying coefficients of v and its derivative in the boundary integral 
constitute the natural boundary condition. That is, 
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dF dF 

t — n x +- — n x =q on T 

au x ally 

is the natural boundary condition. Specification of the dependent variable in the same 
form as the arbitrary test function constitutes the essential boundary condition. In the 
case presented above, only v appears in the boundary integral. Hence, specifying u on T 
is the essential boundary condition. The variables involved in the essential boundary 
conditions of the problem are identified as primary variables and those in the natural 
boundary conditions as the secondary variables in the formulation. The primary variables 
are required to be continuous, whereas the secondary variables may be discontinuous in a 
problem. 

The differential equation is said to describe a scalar-field problem if the 
dependent variable is a scalar and requires only the specification of magnitude for a 
complete description. A vector-field problem is one that requires the specification of 
magnitude and direction. The Poisson equation is an example of a differential equation 
describing a scalar-field problem that arises in many fields of engineering science such as 
elasticity, heat transfer, fluid mechanics, and electrostatics. The equation of motion is an 
example of a differential equation describing the vector-field problem that governs the 
motion of general continua. Each of these categories of differential equations will be 
discussed in more detail and the concomitant formulations presented in Chapters II and 
III. 

The basic concepts of continuum mechanics and the ancillary fundamental 
concepts of mathematical approximation methods outlined in this section form the basis 
for the methodologies developed in this work. In subsequent chapters, the concepts are 
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described further as they relate to the development of multifunctional approaches for 
scalar-field and vector-field problems in engineering science. 

1.3. LITERATURE REVIEW FOR COLLABORATIVE METHODS 

This section includes a literature review of topics related to collaborative methods 
for multi-fidelity modeling and analysis. Review of approaches for collaborative 
modeling of multiple domains is presented. This review is not intended to be an 
exhaustive review of the subject matter but rather to provide sufficient background of the 
fundamental concepts applicable to collaborative methods for engineering science. For 
more detailed discussions on any of the topics reviewed, the reader is directed to the 
referenced reports. 

Multi-fidelity modeling, as referred to herein, entails the use of diverse 
approximations among multiple domains. Numerous approaches for multi-fidelity 
modeling have been developed over the last several decades. Many of these approaches 
are commonplace in the analysis and design of aerospace structures. Generally, these 
methods focus on modeling to obtain accurate stress data, and they have been used 
primarily in an analysis framework rather than as an integral part of the design process. 
With the development of rapid equation solvers and fast computer systems with 
enormous storage capacities, these methods have the potential for impacting the 
preliminary design stage. Research directly applicable to multi-fidelity modeling based 
upon the finite element method continues to flourish. Developments pertinent to this 
research include substructuring, global/local methods, model synthesis methods (i.e., 
multiple method approaches), submodeling, and finite element interface methods. While 
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all of these methods can be used in a global/local analysis, in general, they provide a 
diverse capability for modeling multiple subdomains. 

Substructuring, submodeling, and general global/local methods have been 
highlighted, for example, by Ransom and Ransom and Knight and have been further 
elaborated on by Rose 10 . One notable application of substructuring related to recent 
advances in computational strategies is the use of neural networks to synthesize or 
combine substructures 11 . In reference 11, substructures are modeled individually with 
computational neural networks, and the response of the assembled substructure is 
predicted by synthesizing the neural networks. Statically determinant substructures and 
statically indeterminate substructures were assembled using a supeiposition approach and 
a displacement collocation approach. Typically, substructuring and submodeling 
approaches either require that the finite element nodes along the interdomain boundaries 
coincide or make use of restrictive interpolations of displacements to the boundaries of 
the local models. The global/local method proposed in reference 8 alleviates the 
requirement for nodal compatibility along the local model boundary by introducing a 
surface spline interpolation of the displacements from an independent global model to the 
boundary of a more refined local model. This uncoupled approach was further extended 
to provide global/local model interaction in an iterative approach proposed by Whitcomb 
et al. ’ ' In addition, global/local methodology for two- and three-dimensional stress 
analysis of composite structures has been developed within a common framework by 
Knight et al. 14 

In the context of this work, model synthesis refers to collaborative methodology 
that couples or synthesizes two or more dissimilar mathematical models of multiple 
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subdomains. Myriad methods fall into this category. Examples of these methods and 
representative references include, but are not limited to, synthesis of finite element and 
boundary element methods ’ ’ , finite element and Rayleigh-Ritz approximations , 
finite element and finite difference methods 9 , finite element and analytical solutions , 
and finite element and equivalent plate solutions . Furthermore, an extensive review of 
coupling the finite element method and boundary solution procedures has been given by 
Zienkiewicz . In reference 23, the finite element method is generalized to encompass 
both the finite difference and the finite volume approaches. 

A new era of multi-fidelity modeling was introduced through the development of 
an alternative approach for combining finite element models with different levels of 
fidelity, which is referred to in the literature as interface technology. The concept of 
interface technology is the genesis for the multifunctional capability presented in this 
work. As such, a more extensive review of the literature is presented and the notable 
contributions are outlined. The basic concept of the interface technology was discussed 
by Housner and Aminpour . In this work, the fundamental approaches were discussed 
for mathematically coupling multiple subdomains whose grid points along common 

Of 

boundaries did not coincide. Subsequent developments performed by Aminpour et al. 
implemented the basic concepts, extended the work to alternative approximations, and 
compared the results for representative benchmark applications. Ransom et al. 26 
advanced further the technology by recasting the interface technology in the form of an 
element, thus facilitating the use of the method for more than two subdomains. 

Moreover, the implementation of the method as an element facilitated the inclusion of the 
technology into standard commercially available finite element software codes . Davila 
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et al. 28 extended further the capability for coupling not only along finite element edges as 
originally implemented but across finite element faces as well. Rose 10 extended the 
concept of interface technology to include geometric incompatibility as well as nodal 
incompatibility. In this work, the geometry of the subdomains is automatically adjusted 
to account for an inaccurate geometry description along the common subdomain 
boundaries and for gaps in the boundary definition, which allows for enhanced modeling 
flexibility. In addition, extensions have been developed to include geometrically 
nonlinear analysis . The technology has been developed to provide dimensionality 
reduction for integrating three-dimensional finite element models within two-dimensional 
finite element models 26 . All of the aforementioned interface technology developments 
have focused on a one-dimensional interface along a curve or line. Aminpour et al. and 
Schiermeier et al. have extended the work to a two-dimensional surface interface for 
coupling three-dimensional finite element models. 

1.4, OBJECTIVES AND SCOPE 

The overall objective of this research is to formulate multifunctional methodologies and 
analysis procedures for interfacing diverse domain idealizations including multi-fidelity 
modeling methods and multi-discipline analysis methods. Specific goals of this research 
include: 

1 . To formulate general methodology providing capability for multifunctional 
modeling, analysis, and solution. 

2. To identify computational aspects and related algorithms for this methodology. 

3. To apply the formulation to scalar- and vector-field applications in engineering 


science. 
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The scope of the present work includes the multi-fidelity modeling and analysis of 
interfaced domains within the same discipline as well as among multiple disciplines. The 
analysis capabilities are limited to scalar- and vector-field problems using both single and 
multiple approximation methods within a given domain. The capabilities are developed 
considering discrete changes in domain characteristics across the interfaced boundaries, 
compatibility with general-purpose finite element codes, applicability for a wide range of 
discretization methods and engineering disciplines, and cost-effectiveness related to both 
modeling and analysis time. To accomplish the objectives of the present work, numerical 
studies are performed to gain insight into the interactions among the interfaced domains 
and the computational strategies for the modeling and analysis. Prior to applying the 
method to vector-field problems, the proposed method is evaluated with regard to 
accuracy and computational implications on representative scalar-field problems. 

The organization of the remainder of the dissertation is as follows. A 
multifunctional approach for scalar-field problems is presented in Chapter II. Single- and 
multiple-domain formulations are presented in the chapter along with a discussion of the 
spatial modeling and the computational implications, and numerical results for a 
verification test case are presented. The multifunctional approach for vector-field 
problems is presented in Chapter III. Single- and multiple-domain formulations are 
presented in this chapter along with a discussion of the spatial modeling and the 
computational implications, and numerical results for a verification test case are 
presented. Numerical results for representative scalar-field problems in engineering 
science are presented in Chapter IV, while results for vector-field problems are presented 
in Chapter V. In addition, a discussion of extensions of the methodology to multiple 



discipline coupling is given in Chapter V. Conclusions and recommendations are 
presented in Chapter VI. An overview of the steps in analysis and simulation is given 
Appendix A. A derivation of the cubic spline interpolation matrices used in the 
multifunctional approach is presented in Appendix B. Details of the geometry 
representation along the subdomain interface are given in Appendix C. 



21 


CHAPTER II 

MULTIFUNCTIONAL APPROACH FOR SCALAR-FIELD 

PROBLEMS 

2.1. GENERAL 

The motivation for the consideration of multifunctional approaches for scalar- 
field problems comes from the fact that methods of approximation such as Ritz, Galerkin, 
and other weighted residual methods are based on weak statements of the differential 
equations governing the system response. The differential equation is said to describe a 
scalar-field problem if the dependent variable is a scalar and requires only the 
specification of magnitude for a complete description. The scalar-field problem is a basic 
form of the governing differential equations and thus lends itself to forming the 
mathematical foundation for the general methodology developed herein. Representative 
examples of the scalar-field differential equations in two dimensions are considered 
herein, and the mathematical statement is formulated. The concepts developed here are 
directly applicable to one-dimensional scalar-field problems; however, the development 
is not included in the interest of brevity. The general form of the differential equation 
describing a scalar-field problem for domain Q (see Figure 2.1) is given by the Poisson 
equation, which is of the form 

-V-(XVu) = Q inQ (2.1) 

subject to the natural boundary condition, k— + h(u -u oa ) = q on T , and essential 

dn 

boundary condition, u = u on F P The normal derivative, — = ——n x + — n v , and n x 

dn ox ay 7 
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and n y are the components of the outward normal vector, n, to the bounding surface, T, of 
domain, Q. In Eq. (2.1), the variables k and Q are known coefficients, and the primary 
variable or dependent variable is u, which is a function of the independent variables x and 
y. In the natural boundary condition, the variables, h and u^, are the convection 
coefficient, and the far-field value of the primary variable, respectively. The terms, q, 

k — , and k— are the secondary variables that may be described on a portion of the 
ax ay 

boundary, T . The primary variable, u, is specified on the boundary, T P , and its 
prescription to the boundary value, u , constitutes the essential boundary condition. The 

complete boundary is defined as T = + T s . 



Figure 2. 1 . Geometric Representation of Two-Dimensional Domain. 

2.2. DISCIPLINE SPECIFICS 

Equations of the type of Eq. (2.1) arise in many fields of engineering science such 
as elasticity, heat transfer, fluid mechanics, and electrostatics. Reddy has tabulated 
several examples. In this work, the Poisson equation is applied to problems in the solid 
mechanics and fluid mechanics disciplines. 
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2.2.1. Solid Mechanics 

For applicability of the Eq. (2.1) in solid mechanics, consider a prismatic bar of 
constant cross section subjected to equal and opposite twisting moments at the ends as 
shown in Figure 2.2(a). 



(a) Geometry 


(b) Partial End View 


Figure 2.2. Geometric Configuration of Prismatic Bar. 


In general, the cross sections normal to the axis of the bar warp. As a fundamental 
assumption, the warping deformation is taken to be independent of the axial location and 
is given by 

w= w{x,y) 

Assuming that that no rotation occurs at the endz=0 and that the angle of rotation, 0, is 
small, the displacement components, u and v, in the x andy coordinate directions, of an 
arbitrary point, P, P(x, y), in a plane for constant z, are respectively, 

u = -(rOz)sincc = -yOz 


v = {r6z) cos a = xOz 


( 2 . 2 ) 
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where the angular displacement of a line segment, OP, from the origin, O, to an arbitrary 
point, P, is 0z and a is the angle between OP and the x axis (See Figure 2.2(b)). By 
substituting Eq. (2.2) into the strain-displacement relations, the following are obtained 


'Fv Yxy t'y 0 


dw . dw . 

Y zx — y @ > T^zy — 0 + X® ( 

The three-dimensional stress-strain relations given in terms of Lame’s constants for a 
linear isotropic solid are given by 

(7 x — 2Ge x + Yc , Try — Gy xy 
@ ' y 2Ge y 4" %£■ , x ’ y2 Gy y2 
o z — ^-Ge z + A& , t xz — Gy xz 


where e - e T + E v + e 7 , X = w r , and G = r . 

x 7 z (l + v)(l - 2v) 2(l + v) 

The shear modulus, G, and the quantity, X, are referred to as the Lame’s constants, and 

the modulus of elasticity, E, and the Poisson’s ratio, v, are material properties. 

Substituting the strain-displacement relations of Eqs. (2.3) into the stress-strain relations 


— T xy ~ ~ ~ 0 


r„ =g(— -yfll ; t„=gI \—+ x6 

la* 7 * Uv 


Then, the three-dimensional equations of equilibrium, 
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+ - 


d(J Y dr 
dx 
dx 


xy 


dx 

^ Tvz 


+ - 


dy 

dOy 

dy 

dx 


x^ + dT 


- + 


+ F y = 0 

dz x 

dX yg 

— — + F v = 0 

dz y 


(2.5) 


+ - 


dx dy 


yz d(J 
' + 


dz +Fz 0 


with negligible body forces, simplify to the following equations: 

dx dx 7V r)x dx?,, 

^SL = 0, -21 = 0, ^+—^ = 0 (2.6) 

dz dz dx dy 

First, note that the stresses in Eq. (2.4) satisfy exactly the first two equilibrium equations 
above (see Eq. (2.5)). Next, Eq. (2.4) can be combined into a single equation by 
differentiating the expressions for x zx and x zy by y and x, respectively, and subtracting the 
resulting equations. These operations yield the compatibility equation given by 

dx zx 


= - 2 G0 . 


(2.7) 


dy dx 

The stress in a bar of arbitrary cross section may thus be determined by solving the third 
equation of equilibrium given in Eq. (2.6) along with the equations of compatibility given 
in Eq. (2.7) and the given boundary conditions. 

This torsion problem is commonly solved by introducing a single stress function. 
If such a function, <|)(x, y ), the so-called Prandtl stress function, is assumed to exist, such 
that 


r =<!£ • t =-$£ 

zx dy ’ ^ dx’ 

then, the equations of equilibrium are automatically satisfied. The equation of 
compatibility becomes, upon substituting these expressions for the shear stress. 
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d 2 (f) | d 2 (t) 
dx 2 + dx 2 


= -2 GO 


Therefore, if the compatibility requirement is to be satisfied, the stress function, (|), must 
satisfy Poisson’s equation, Eq. (2.1). The primary variable, u, the constant, k, and the 
source variable, Q, are represented by the stress function, the inverse of the shear 
modulus, G, and twice the angle of twist per unit length, 0, respectively. Moreover, the 
stress function, (|)=constant on the surface of the bar. 

2.2.2. Fluid Mechanics 

For a two-dimensional incompressible irrotational flow, expressions are given for 
the velocity components, v x and v y , in terms of the x andy coordinate directions, 
respectively. The velocity components should satisfy the continuity condition 


dv x OVy 

V • v = — — + — — 

dx dy 


and the irrotational flow condition 


dx dy 


In terms of the stream function, \|/, the components are given by 


dy , ay 

v,=— and v y = — 


and in terms of the velocity potential, O, the components are 


ao a ao 

v r = and v v = 

x dx y dy 


Substituting the velocity components, v x and v y , from Eq. (2.10) into the irrotational flow 
condition Eq. (2.9), one obtains 
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d 2y V d 2x ¥ A 

^ + r- = 0 . 

2 'v . 2 


( 2 . 12 ) 


dx z dy ' 

Note that the velocity components in terms of the stream function given in Eq. (2. 10) 
satisfy the continuity condition, Eq. (2.8) identically. Hence, Eq. (2. 12) governs the flow 
in terms of the stream function, \|/, and is in the form of the Poisson Equation, Eq. (2.1) 
where the primary variable, u, the constant, k, and the source variable, Q, are represented 
by the stream function, \|/, the density, p, and the mass production, a (normally zero), 
respectively. 

Substituting the velocity components, v x and v y , from Eq. (2. 1 1) into the 
continuity equation, one obtains 


a 2 o a 2 o . 

-t + tt =0 - 


(2.13) 


dx z a/ 

Note that the velocity components in terms of the velocity potential given in Eq. (2.1 1) 
satisfy the irrotational flow condition, Eq. (2.9), identically. Eq. (2.13) governs the flow 
in terms of the velocity potential, <t>, and is in the form of the Poisson Equation, Eq. (2.1), 
where the primary variable, u, the material constant, k, and the source variable, 0. are 
represented by the velocity potential, O, the density, p, and the mass production, o 
(normally zero), respectively. 

2.3. SINGLE-DOMAIN FORMULATION 


In this section, multifunctional methodology for a scalar- field problem over a 
single domain is presented in terms of weighted residuals. The method of weighted 
residuals is used extensively in fluid mechanics and thus the potential problem is 
formulated from this perspective. While the intent of this work is to develop general 



28 


methodology for multiple domains, the salient features of the weighted residual method 
formulation may be investigated and discussed using the single domain. Consider the 
general Poisson equation for a two-dimensional domain for field variable, u 

-kV 2 u = Q (2.14) 

in a domain, Q, bounded by T. In general, the boundary, T, can have mixed boundary 
conditions with the primary variable, u, prescribed on F and the secondary variable, the 

S 

flux, q, prescribed on the remaining part of the boundary, T (see Figure 2.1). 

In the method of weighted residuals, an approximate solution, u , is used in 

expressing V u , then the differential equation, Eq. (2.14), will no longer be satisfied, 
and this lack of equality is a measure of the departure of u from the exact solution. The 
lack of equality is called the residual, R, and is written as 

R = -kV 2 u — Q^ 0 . 

The residual is orthogonalized by a set of weight functions, O, and averaged over the 
domain. This residual may be written as 

l(-kV 2 u-Q) O / dn = 0. (2.15) 

Q 

n 

The solution for u is sought in the form u = Y. a i^i + T'q . The functions, *F,, are usually 

i=l 

called trial functions, and a,- are arbitrary constants. The trial functions satisfy the 
homogeneous boundary conditions, while 'To satisfies the nonhomogeneous boundary 
conditions. Posing the problem to be solved in a generalized weighted residual form 33,34 
and relaxing the requirement for the approximate solution to satisfy all boundary 
conditions, the weighted residual statement may be written in the form 
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J<PA(5)dQ + <f<I>B(5)dr = 0 

o r 

where the residual in the satisfaction of the boundary conditions is orthogonalized by a 
secondary set of weight functions, <t> , and the differential equation set is represented by 

A x {u) 

A(u) = -i ^(w) ► = 0 

in the domain, Q, together with the boundary conditions 

B l(«) 

B (“) = 1 5 2(w)r = 0 
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Note that the trial function may be selected so as to satisfy the essential and the natural 
boundary conditions; thus, the boundary integrals in Eq. (2. 16) are identically zero. In 
this formulation, only the essential boundary conditions, i.e., 


u - u = 0 on r p 

are assumed to be automatically satisfied by the choice of the trial functions. Therefore, 
Eq. (2. 1 6) is rewritten as 



(2.17) 


or 



r d 2 u d 2 u ' 



r d u _) 

-k 


-Q 

dD+ ftp 

k- q 


dx 2 dy 2 
\ 7 ) 


J 

r s 

l dn J 


Idr 5 = o 


H 

Q 


. du du du , ^ — 

where — = — n Y H n v and <Po = ‘I> . 

d n dx x dy y 2 

In general, the method of weighted residuals does not strictly require the 
incorporation of natural boundary conditions into the weak formulation, as in the Ritz 
method. However, if the operator permits the weak formulation, continuity requirements 
on the primary variable and its derivatives may be relaxed. Moreover, if integration by 
parts is possible, one may reduce the order of the highest derivative in the integral form 
to eliminate the difficulty of selecting the appropriate weight functions. Thus, in the 
formulation herein, the order of differentiation on the primary variable in the integral 
equation, Eq. (2. 17), is reduced to obtain the weak formulation. In addition, 
acknowledging that the primary variable, u, is approximated by u , for simplicity, the 


subsequent development is presented in terms of u. Application of the divergence 
theorem to Eq. (2.17) yields 
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\ k \ 

£2 


du 90 du 90 
dx dx dy dy 


■ r , | du du 

bo ' H ¥"' + ^ 


bar 


(2.18) 

-JgOaT2 + \ \k— -q fear 5 =0 
£2 pA ^ ) 

Note that the bounaary is presumea to consist of bounaaries on which the primary 
variable is specifiea ana bounaaries on which the seconaary variable is specifiea, ana 

r = F^ +T S . Therefore, the bounaary integral on V. may be expressea as 


du du 


dy 


r du du ^ 


jk\ f x n x + %n y |oar= J/c| + %n y fMT* + J *| 


T p 


dx dy 


Ts 


du du Y_ 
—n Y +—n v bar 


9x 


dy 


Noting that, in the methoa of weightea resiauals, the weight function, O, satisfies the 
homogeneous bounaary conaitions for the primary variable (i.e., essential bounaary 
conaitions). Thus, 0=0 on r . Therefore, the bounaary integral on l 4 is iaentically zero 
ana Eq. (2. 1 8) may be rewritten as 


l /f l 

£2 


du 30 du 30 
+ 




_ . ^ . .. oar 5 - /0oaD+ <f 

dx dx dy dy y pS dn Q 


, du _ 
k q 

dn 


b dr 5 = 0 


Since the weight functions, O ana O , are arbitrary, they may be chosen, without loss of 
generality, such that, 0 = 0. Therefore, 


\ k \ 

£2 


du 30 du 30 
dx dx dy dy J~ q 


j£>oaQ- f^oar 5 = 0 


or 


l /c l 

£2 


du 30 du 30 


• + - 


dx dx dy dy 


= jQ OdQ+ jq® ar ? 


£2 


(2.19) 


The integral form of Eq. (2.19) forms the basis of finite element approximations, which is 
summarizea in a subsequent section. 
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2.4. MULTIPLE-DOMAIN FORMULATION 

In the multiple domain method, the domain of the problem is subdivided into a 
number of smaller subdomains. The method is quite similar to the subdomain collocation 
method, which is another weighted residual method. In the subdomain collocation 
approach, the domain is divided into as many subdomains as there are adjustable 
parameters. These parameters are then determined by making the residual orthogonal to 
a weight function in an integral sense over each subdomain. Here, as in the single- 
domain formulation, methodology is presented formulating the general method of 
weighted residuals for multiple domains by considering the Poisson equation for a two- 
dimensional domain for a field variable, it. Then, 

-kV 2 u = Q (2.20) 

in the entire domain, Q, bounded by 1 . For simplicity, the multiple-domain formulation 
is presented for two subdomains, Qi and Q 2 (see Figure 2.3). Independent 
approximations and weight functions are assumed in each of the subdomains and 
continuity conditions are used to provide for a continuous solution across the subdomain 
interfaces. Thus, Eq. (2.20) is satisfied in each subdomain, independently, i.e., 

-k{V 2 U\ =Q\ in and -kfl 2 Ui = Q 2 i n ^2 
subject to the boundary conditions on the subdomain boundaries, Ti and T 2 . Although 
Eq. (2.20) is assumed for uniform constant, k, throughout the domain, it is permitted to be 
different in each subdomain. That is, constants, k\ and fe, are used for subdomains Qi 
and Q 2 , respectively, to allow for the general case of nonhomogeneous material. 
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At this point, differences between the single- and multiple-domain approaches 
become evident. First, the domain, Q, is now represented by the union of ns subdomains, 
Q/, such that 

ns 

Q. = . 

i'=l 

Second, the bounding surface, T, of the domain, Q, is the union of the exterior surfaces, 
r, , of the ns subdomains, Q t , such that 


ns 


r=srf 


i = 1 


p 

In general, these exterior surfaces, p , may involve mixed boundary conditions with the 
primary variable, u, prescribed on r f and the secondary variable, the flux, q, prescribed 
on r f such that 


rf =rP +Tf . 

Finally, as a result of the subdomain modeling, the collaborative effort to solve the 
problem involves an interior surface interface boundary, r/ , and the information transfer 
across the boundary. Hence, the boundary surface for the i th subdomain is given by 

r ( =r ( p +rf + r/ 

The boundary conditions may be written as 

U\ - U\ = 0 on r/ and k\ — - -q\=0 on Ff 

d n 


and 
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U 2 ~U 2 =0 on 


rP 

l 2 


and 


d u 2 - 
k 2 ~ q, 

an 


= 0 on 


y s 

l 2' 



Figure 2.3. Boundary Definitions for Two-Dimensional Subdomains. 


The residual for each domain is orthogonalized by a set of weight functions, <E>, and is 
written as 


j(-*lV 2 «i -0i )oj dD 1= 0 


and 


j(-^ 2 V 2 W2 -£>2)^2 - 0 

Q 2 


^ n 

where the approximate solution is sought in the form u\ = Y, a \i^\i + ^0/ an d 

1 


^ n 

U 2 = Y, a 2i^2i + VQi * The functions, T 0 / and ¥ 2 /, are trial functions, and au and c/ 2 / 
1 
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are sets of arbitrary constants. Using the general form outlined in the single-domain 

formulation (/. e. , JOA(u)dQ + f<PB(u) dr = 0 ), for each subdomain, one may write 
o r 

JO ; A(S ( -)dO+ «fO f B(5 f )dr, =0 for i = 1,2 
U 

where the residual in the satisfaction of the boundary conditions is orthogonalized by a 
secondary set of weight functions, O, , for subdomain i. Therefore, considering the 
approximate solution, u\ and > we may write the general integral form of the 
differential equation governing the potential flow for subdomain 1 as 

JO 1 (-£ 1 V 2 M 1 -0 1 )dQ 1 + V* =0(2.21) 

Qi r p r° V m J 

and for subdomain 2 as 


J ® 2 ^ 2 ^ ^2 ~ Q 2)^2 + \^ > 2 \& 2 ' u 2) ( ^2 + J < ^ > 22 ^ 2 ' 

rf rf v 


-q 2 dr" = 0 ( 2 . 22 ) 


Again, we will presume that the essential boundary conditions, i.e., 


and 


u\ - U\ = 0 on Tf 


«2 ~ u 2 = 0 on ^2 

are automatically satisfied by the choice of the functions, U\ and . Therefore, for 
subdomain 1, Eq. (2.21) is rewritten as 



(2.23) 


or in its expanded form 
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-k 


f d% d%} 

dx 2 + dy 2 


-Qi 


dDi + J <I>] 


r / 


d u\ _ 

k \~ A 9\ 

an 


dr/ = o 


(2.24) 


where— — = —^n x . +^- L n y . and <I>i 2 = ®i and ®22 = ®2- Similarly, the weighted 
dn dx 1 dy 


dui 


residual form for subdomain 2, 


J®2 

q 2 


( 2 ~ 2 ~ 
o U2 a U2 

dx 2 + dy 2 


\ 

_ ( 

-02 

dD 2 + J ®2 

J 

r 2 ^ 


dw? _ 


dn 


■92 


dTf = 0 (2.25) 


The order of differentiation on the primary variable in the integral equations, Eq. 
(2.24) and (2.25), is reduced to obtain the weak formulation. In addition, acknowledging 
that the primary variables, U\ and 7 / 2 , are approximated by U\ and «2 » f° r simplicity, the 
subsequent development is presented in terms of U\ and 7 / 2 . Utilizing the divergence 
theorem, Eq. (2.24), can be rewritten, for subdomain 1, yields. 


J h 


y du\ dOj du\ d®| t , ( dwi d u 


■ + ■ 


dx dx dy dy 


n 


\ 


1 OMl 

— L n Y , h L n v , 

v dx x > dy yi y 


Ml- Jft^l 


+ 


d u\ _ 

h-r~9\ 

d/? 


(2.26) 


|<I>ldrf = 0 


and similarly, for subdomain 2, 
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(2.27) 

Note that the domain boundary is presumed to consist of boundaries on which the 
primary variable is specified, boundaries on which the secondary variable is specified. 
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and boundaries at the subdomain interface. Thus, for subdomain i, r, = Yf + Yf + Yj . 
Here, the boundary on the interface is assumed to be conforming (i.e., represents same 


geometry) and F* = f 1 . Therefore, the boundary integral on F/, may be expressed as 


Ff 


duj V f ( 

— " “ P/dTf=J 
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Note that <I>, = Oon T, . Eq. (2.27) can be rewritten, for subdomain 1, as 
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Since the weight functions, ®] and ®] , are arbitrary, they may be chosen, such that, 
<I>i = d>| . Therefore, 


J h 


o, 


du\ d<!>i dii\ 


ar ar - ^ 

J T T { 

- Jg^dD! = 0 
«1 


(2.28) 


Similarly, for subdomain, Q 2 , 
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(2.29) 
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In general, for the multiple domain case, the approximation for the primary 
variable ( e.g the potential field) must satisfy the following conditions: 

i. The primary variable must be continuous and single valued in the subdomain. 

ii. The primary variable must be continuous across the interdomain boundary. 

Hi. The primary variable on the subdomain boundary must satisfy the boundary 

conditions. 

If the requirement to satisfy interdomain continuity is relaxed, an additional boundary 
condition is used, namely, 


U\ — u 2 = 0 on r . 

This constraint can be satisfied in the integral sense as 

JA(«! - w 2 )dr = 0 on r 1 (2.30) 

Ti 

where X is a Lagrange multiplier associated with the secondary variable along the 
common subdomain boundary. Therefore, combining Eqs. (2.28) and (2.29) for the 
entire domain, and including the continuity integral at the interdomain boundary yields 
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where for subdomain, i, q,- , are the secondary variables along the interdomain 


, , » d Uj dui duj 

boundary, q, = — - = — L n x + — L n 
on ox 1 dy 


yt 


Note that the normals on the inter domain 
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boundary are equal and opposite (see Figure 2.3). That is, n j = -n 2 where 
n ; = n x . i + fly. j , and it follows that, qj = -q 2 = q . Therefore, 
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or rearranging 


J k \ 


Qi 


du\ 3<1>[ du\ 30[ 


dx dx dy dy 


IdQj - jq < Dj dr 1 - JftOjdQ! - \q x ^\ dTf 


Qi 


+ 


J k 2 

On 


r du 2 3 ® 2 + du 2 d<S>2 V, , is*, jrl 


^ . .. . m 2 + dT - }02<M^2 - J?2®2 dP 

d* dx dy dy J p i q 2 r * 


+ 


J/l(iq -w 2 )dr : 

r^I 


= 0 


(2.32) 


Note that Eq. (2.32) is written as a single equation for convenience and represents the 
sum of terms related to the residual in the governing differential equation within each 
subdomain and the continuity constraint for the primary variables along the common 
subdomain boundary. However, each of the bracketed terms in Eq. (2.32) must equal 
zero individually. These bracketed terms are identical to Eqs. (2.28), (2.29), and (2.30) 
which must be satisfied independently. 

In this formulation, the two primary field variables, U\ and u 2 are approximated 
independently, and continuity requirements between these two approximations are 
satisfied along the subdomain interface boundary. The use of these approximations and 
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the associated continuity requirements gives rise to the classification of the formulation 
as a two-approximation approach. 

Now consider a configuration that makes use of a third approximation for the 
primary variable along the subdomain interface boundary in addition to the 
approximations given along the boundary of each of the subdomains. This primary 
variable, v, along the interface is assumed to be independent of the primary variables, U\ 
and « 2 , of the subdomains to which it is attached. These independent approximations 
give rise to continuity requirements along the interface of the form 

v — u\ = 0 on r 1 

v-u 2 = 0 on r 1 

These constraints can be satisfied in the integral sense as 

/^(v-M^dr 1 =0 on r 1 (2.33) 

r 1 

\X 2 (v-u 2 )&T l =0 on T 1 (2.34) 

r 1 

where X\ and X 2 are Lagrange multipliers or weight functions in the foim of the 
secondary variable along the interface. An additional continuity requirement in terms of 
the secondary variable along the common subdomain boundary is required. These 
secondary variables, q\ and q 2 , are assumed to be independent of each other. These 
independent approximations give rise to continuity requirements along the interface of 
the form 

q\+q 2 =$ on T 1 
These constraints can be satisfied in the integral sense as 
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\X{q x +q 2 )dT l =0 on T 1 

T^I 


(2.35) 


where X is a Lagrange multiplier or weight function of the form of the primary variable 
along the interface. Combining Eqs. (2.28) and (2.29) for the entire domain and 
including the three continuity integrals along the interdomain boundary, Eqs. (2.33), 
(2.34), and (2.35), yields 
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Again, note that Eq. (2.36) is written as a single equation for convenience and represents 
the sum of terms related to the residual in the governing differential equation within each 
subdomain and the continuity constraints for the primary and secondary variables along 
the common subdomain boundary. Each of the bracketed terms in Eq. (2.36) must equal 
zero individually. These bracketed terms are identical to Eqs. (2.28), (2.29), (2.33), 
(2.34), and (2.35), which must be satisfied independently. 
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The integral form of Eqs. (2.32) and (2.36) forms the basis for the subsequent 
spatial modeling approximations. The spatial modeling approximations are discussed in 
detail in the next section. Eqs. (2.32) and (2.36) may be generalized for more than two 
subdomains and for multiple interfaces by 
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(2.38) 


where N ss is the number of subdomains in which the entire domain is subdivided, Nj is 
the number of interfaces connecting the N ss subdomains and n ss (i) are the number of 

subdomains attached to interface i. For example, for one interface connecting two 
subdomains, Eq. (2.38) yields in its expanded form 
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which is identical to Eq. (2.36). 


2,5. SPATIAL MODELING FOR MULTIPLE DOMAINS 

Although this section is focused on spatial modeling of multiple domains using a 
multifunctional development, a brief discussion of spatial modeling for a single domain is 
given first, followed by a more detailed discussion for multiple domains. Thus far, a 
multifunctional approach based on weighted residuals has been formulated. This 
approximation technique provides a mechanism for finding approximate solutions to 
problems in mathematical physics and engineering science such as those represented by 
the Poisson problem. Selection of the approximating and weighting functions for 
complex geometrical shapes and boundary conditions poses a major difficulty for 
weighted residual methods. In addition, the methods were generally not regarded as 
being computationally competitive compared to the traditional finite difference method. 
However, weighted residual methods offer a versatile means by which to formulate finite 
element equations where no functional is available. Hence, many of the difficulties 
associated with this class of methods are alleviated. The derivation of discrete equations 
is an essential component of the approximation technique. Thus, several discretization 
approaches are outlined in the next section. 
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2.5.1. Overview of Discretization Methods 

Various forms of spatial modeling or discretization of the continuum problem 
defined by the differential equations can be used. These forms include, but are not 
limited to, the finite difference method, the finite volume method, the finite element 
method, and the boundary element method. In such spatial modeling, the infinite set of 
numbers, representing the unknown function or functions is replaced by a finite number 
of unknown parameters. A brief discussion of each of the aforementioned modeling 
methods is given here to provide the foundation for discussion of interfacing such diverse 
methods, which is presented in subsequent subsections. 

The finite difference method 

Of the various forms of spatial modeling, one of the simplest is the finite 
difference method. The finite difference method gives a pointwise approximation to the 
governing equations. In the finite difference approximation of a differential equation, the 
derivatives in the equation are replaced by differential quotients that involve the values of 
the solution at discrete mesh points of the domain. The resulting discrete equations are 
solved for values of the solution at the mesh points, after imposing the boundary 
conditions. While finite difference techniques are widely used in fluid dynamics and heat 
transfer and can treat fairly difficult problems, they become hard to use when irregular 
geometrical shapes or unusual boundary conditions are encountered. In addition, because 
it is difficult to vary the size of different cells in particular regions, the method is not 
suitable for problems of rapidly changing variables, such as stress concentration 
problems. These adverse attributes are particularly significant in structural analysis. 
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The finite volume method 

The finite volume method evolved in the early seventies via the finite difference 
approximations and has many proponents in the field of fluid mechanics. The method 
takes as its starting point the physical conservation laws in integral form written for small 
control volumes around every discrete point. Modifying the shape and location of the 
control volumes associated with a given discrete point, as well as varying the rules and 
accuracy for the evaluation of the fluxes through the control volume, gives the method 
considerable flexibility. Unlike the finite difference method, the finite volume method 
can readily handle arbitrary mesh orientation thus making it more amenable to problems 
of rapidly changing variables. In addition, by direct discretization of the integral form of 
the conservation laws, the basic quantities (e.g., mass, momentum, and energy) will be 
conserved at the discrete level. Like the finite difference method, the finite volume 
method has been shown to be a special case of the finite element method with non- 
Galerkin weight functions . 

The finite element method 

The finite element method consists of representing a given domain by an 
assembly of smaller, geometrically simple subdomains or elements over which the 
approximation functions are systematically derived. Then, Ritz-Galerkin approximations 
of the governing equations are developed over each element. Finally, the equations over 
all elements of the collection are connected by continuity of the primary variables. In the 
mathematical literature, the names Petrov-Galerkin are often associated with the use of 
weighting functions such that O ^ N , and the names Bubnov-Galerkin are often 
associated with the use of weighting functions such that <1> = N , where in the finite 
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element method N are the element shape functions. The latter method is often referred to 
as the Galerkin method. The resulting system of equations is sparse, banded, symmetric, 
and positive definite. The finite element method is especially well suited for handling 
arbitrary shapes or domains. To obtain good accuracy in regions of rapidly changing 
variables a large number of small elements must be used. Furthermore, the method is 
widely used for the analysis of many engineering problems involving static, dynamic, and 
thermal stresses of structures. 

The boundary element method 

The boundary element method is an alternative to the finite element method. 

Like the finite element method, the boundary element method uses nodes and elements to 
discretize the boundary of the domain. Thus, compared to the finite element method, the 
dimensionality is reduced by one. The governing differential equations are transformed 
into integral identities, which are applicable over a surface or boundary. These integrals 
are numerically integrated over the boundary, which is divided into small boundary 
segments. The method may be used to model accurately the response in the domain 
bounded by its mesh. The method can easily accommodate geometrically complex 
boundaries. Furthermore, since all the approximations are restricted to the surface, the 
method can be used to model regions with rapidly changing variables with better 
accuracy than the finite element method. Complex kernel routines are required to 
determine the response for the interior of the domain. Hence, the computational expense 
increases quickly if the response at several interior locations is needed. In addition, for 
nonlinear problems, the interior must be modeled; thus losing the advantage of reduction 
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in dimensionality. Unlike the finite element system matrix, the analogous boundary 
element matrix is small, fully populated, and unsymmetric. 

Each of the aforementioned discretization approaches has advantages and 
disadvantages specific to the domain of the physical problem or the discipline within 
which it is applied. To overcome the disadvantages of the individual methods, coupled or 
collaborative methods have been developed. Collaborative methods couple two or more 
discretization approaches and make use of a given approach when and where it is best 
suited. The interaction between the methods is an essential feature related to the 
robustness and accuracy of the combined methods and is a subject of discussion herein. 
Moreover, this work focuses on the application of the multifunctional method developed 
here to the finite difference and finite element methods and their coupling. 

Computational methods using finite-differences for fluids experiencing field 
discontinuities such as shock-waves and flow separations have been proven to be 
efficient solution techniques. The finite element method has proven to be efficient in 
solving for the response of complex aerospace structures, which may contain internal 
discontinuous members such as spars, ribs, and bulkheads found in fuselage and wing 
structures. In addition, coupled finite difference/finite element methods have been 
proposed that make use of the strengths of the each of the modeling methods in the 
solution of the aeroelastic problem and elasticity problems in references 36 and 19, 
respectively. Thus, both spatial modeling approaches and their coupling will be 


discussed in turn. 
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2.5.2. Overview of Single-Domain Spatial Modeling 
Finite element discretization 

For a single domain, the finite element equations may be obtained by rewriting 
Eq. (2.19) over an element domain as 
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a e 


du d<I> du d<I> 

dx dx dy dy 


= |0<DdD e + jqOdr s 


(2.39) 


Q 


where superscripts on the domain, Q, and boundary surface, T , integrals denote 
integration over the element. In later sections, numeric subscripts will be used to denote 
element integration within the specified subdomain. The primary variable is 
approximated over the element domain by u = Nu g , and using the Galerkin method, the 

vector of weight functions is given by <I> = N . Substituting approximations into the 
integral equation given in Eq. (2.39) yields 
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or 

k e u e =f e 

where k e is the element stiffness matrix, u, is the vector containing the generalized 
primary variables, and i e is the element force vector containing the generalized secondary 
variables. The element field quantities, k, u, and f, are denoted by a subscript, e. 
Assembling these element equations over the entire domain and enforcing continuity of 
the primary variable at the interelement boundaries yields the system of equations given 
by 


Ku = F 



49 


nelem 

where K = Y J k\ 

e=l & 


dN dN T dN 
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dx dx dy dy 


ml ; u is the assembly of all of the nodal 


degrees of freedom associated with the primary variables; and 


nnodes 


F = X jN r 0dQ e + \N J qdT s 
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Finite difference discretization 

In the finite difference methods, derivatives are approximated by difference 
expressions that transform the derivatives appearing in the partial differential equations to 
algebraic equations. For an elliptic partial differential equation, usually time- 
independent, the methods result in a system of algebraic equations that are solved using a 
direct or iterative solution technique. For hyperbolic and parabolic partial differential 
equations, a set of algebraic equations is obtained. These equations are solved either 
explicitly or implicitly. For the explicit solution, each equation will yield one unknown. 
The matrix of unknown variables is a diagonal matrix and the right-hand-side vector of 
the system is dependent on the variables at previous times. For the implicit solution, the 
equations are coupled and must be solved simultaneously. Since the system equations are 
coupled and more than one set of variables is unknown at the same level, the matrix to be 
inverted is non-diagonal. In most cases, however, the structure of the matrix will be 
rather simple, such as a block pentadiagonal, block tridiagonal, or block bidiagonal. The 
truncation errors, stability and consistency of the numerical scheme are aspects that must 
be considered in the development of the methods. The difference expressions are 
obtained by Taylor series expansion, using forward, backward or central expansions. 
Zienkiewicz and Morgan have shown that the finite difference method of approximation 
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is a particular case of collocation with locally defined basis functions. In the collocation 
method, the unknown weight function parameters are determined by forcing the residual 
in the approximation to vanish at N selected points in the domain. Upon substitution of 
the approximation function into the differential equation, the equations can be recast in 
weighted residual form by selecting <I>, = Six -Xj). The weighted form of the residual 

reduces to the evaluation of the partial differential equations using the approximate 
solution evaluated at the N selected mesh points. For a second-order ordinary 
differential equation, the approximate solution, u , may be given as a function of the 
solution at neighboring points (see Figure 2.4) as 

u=u i - l Nf_ l+ u i Nf + u i+ \Nf + ^ 

where Nf are locally defined quadratic basis functions represented by 



Figure 2.4. One-Dimensional Finite Difference Element Configuration. 
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The locally defined quadratic basis functions, Nf , given here in the Cartesian coordinate 
system, may be written in the element natural coordinate system, c, as 

Nf- i=~f(i-£); Nf=( l-f 2 ); JVf + 1 =if( 1 + 4 ) 

where t,-x/h e . Note that -1 < q < 1 . These basis functions are the standard Lagrangian 
shape functions for three-node one-dimensional finite elements. This derivation for one- 
dimensional problems may be extended to two- and three-dimensional problems. The 
derivation is given for two-dimensional problems considering the bi-quadratic shape 
functions for a nine-node two-dimensional finite element. A schematic of the finite 
difference template and the associated finite element are shown in Figure 2.5 where the 
open circles represent grid points in the five-point finite difference template used to 
represent second-order derivatives. 

The shape functions for a nine-node quadrilateral' 14 are given in Table 2.1. For 
example, the shape function at point i,j - 1 is given by 

Similarly, 

N M ,j = )- - £ 2 1 1 - >7 2 )’ 


and 


N Uj ={^~e\-ii 2 \ 
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The standard finite difference representation follows by direct substitution. This 


specialization of the finite difference method as a form of the generalized method of 
weighted residuals forms the basis for its inclusion in this multifunctional derivation. 



Figure 2.5. Two-Dimensional Finite Difference Element Configuration. 




For a single domain, as in the finite element method, the finite difference 
equations may be obtained by interrogating the weighted residual equations over an 
element domain where the element, e, surrounds node i (see Figure 2.5). The 
approximate solution for the primary variable is given by 

M 

u= £ N m u m or u = Nu e 
m = 1 

where M is the number of shape functions over the element, and the weight function, O,, 
is given by the Dirac delta function, S(x - x ( ,y — y t ) = S(x i ,y i ) . Note that the subscript 

i on the weight function is used to denote the subdomain, while the subscript i on the 
coordinate values, x and y, is used to represent the point in the physical domain at which 
the Dirac delta function is evaluated. Therefore, Eq. (2.39) becomes 
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(See BracewelF’ 8 ), the element equation reduces to 


« e =Q{x i ,y i )+q{x i ,y i ) 

For the second derivative difference approximation, the number of shape functions of an 
element, M= 3 and uf = {u(_\ u t u i+ \}. Therefore, as in the finite element method, 
the resulting finite difference equations may be written in matrix form as 
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where k e is the finite difference “element stiffness” matrix, u e is the vector of generalized 
primary variables and f e is the finite difference generalized force vector. Assembling the 
element equations yields 

Ku = F 

where u contains all of the nodal degrees of freedom associated with the primary 
variables, 


Nelem 

K= I -k\ 
1 


3 2 A 

dx 2 


+ - 


X=Xj 

y=yt 


a 2 iv 

dy 2 


X=Xi 

y=yt 


and 

Nnodes 

F= UQ{xi,yi)+q{xi,yi)\. 
l 

While a single spatial modeling approach (i.e., the finite element method or the 
finite difference method) is used for the single domain formulation, subdomain modeling 
permits multiple discretization strategies to be used in a collaborative manner. These 
discretization strategies include homogeneous approaches in which the same 
discretization method is used in each subdomain and heterogeneous approaches in which 
different discretization methods are used amongst the subdomains. Each of these 
discretization strategies is discussed in the following sections. 
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Table 2. 1 . Shape Functions for a Nine-Node Quadrilateral Finite Element. 



2.5.3. Multiple-Domain Modeling - Homogeneous Discretization 

In this context, homogeneous discretization approaches are applicable to multiple 
subdomain discretization. These approaches make use of a single discretization method 
among all subdomains in which the domain is subdivided. Of the many spatial modeling 
approaches, this work will focus on the finite element and the finite difference methods. 
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Two-approximation interface modeling 

For homogeneous domain discretization developed herein, Eq. (2.32) is used to 
provide the mathematical basis. The two-approximation formulation, for both the finite 
element (FE) and finite difference (FD) methods, may be obtained by rewriting Eq. (2.32) 
over an element domain as 


J k \ 
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Note that the integration over the common subdomain boundary, P, is considered only 
for element edges along that boundary. 

The form of Eq. (2.41) for the two (FE and FD) methods differs by the form of 
the element shape functions and the approximation selected for the weight functions, O. 
For the generalized element expansion of subdomain i, the independent approximations 
for the element generalized primary variables, (i.e . , displacements or velocities), interface 
secondary variables (i.e., tractions or fluxes), and the weight functions associated with the 
secondary variables, are, respectively 

u/ =NjU e . ; q = R,« and Xj = R, 


where « is a vector of unknown coefficients associated with secondary variable, q , and 
N and R are matrices of interpolation functions for the element primary and secondary 
variables, respectively. The interpolation functions in the matrix R are assumed to be 
constants for linear finite elements and linear functions for quadratic finite elements. 
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Substituting these approximations into Eq. (2.41) yields an integral equation in terms of 
the weight function, <I> , which is given by 
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f e/ = Jo;adnf+ dr/ e 
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Assembling the element equations over the entire domain, enforcing continuity of the 
primary variable only within each subdomain and assembling the contributions along the 
element edges on the common subdomain boundary, and noting that u C| and u e2 and 


f C] and f e2 are completely uncoupled, yields the system of equations given by 
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0 K, 
Ko K c 


(2.42) 


The system of equations given in Eq. (2.42) is obtained based on the initial development 
of the weighted residual statement, from Eqs. (2.28) and (2.29), 

f , fdUj 9<I>; du,- 3<I>; , , d Us ^ T , f _ _ 

J h 1 k i~r^i dr - + Si&af 

Q . ^dx ox dy dy J r i d n Q . r? 


subject to the constraint equation, Eq. (2.30), 


\X{u x -u 2 )dT l =0 on T 1 . 


Here, the first two matrix equations in Eq. (2.42) are obtained from the weighted residual 
statement for each subdomain, Eqs. (2.28) and (2.29), and the third matrix equation is 
obtained from the constraint on the primary variables along the common subdomain 
boundary, Eq. (2.30). 

For the finite element modeling, the weight functions are taken to be the finite 
element shape functions. That is, <I> ( = N,- , and thus, for i=l,2 


ke , J JwZjHt + wLWl W , 

' L Bx Bx By By 
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Ut= J N/O/dQf + jN^-dr/ 


Q ? if 


Here, note that at the element level, k p = ' , and consequently, at the global system 


level, K„ =K T . 

Pi Sj 


For the finite difference modeling, the weight functions are taken to be the Dirac 
delta function. That is, <X> f = S t (x-x^y- y t ) = S t {xj,yj ) , and thus, for i=l ,2 
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Three-approximation interface modeling 

For the three-approximation formulation, Eq. (2.36) is used to provide the 
mathematical basis for the development. In previous work by Aminpour et al. , a 
similar formulation based on the principle of minimum potential energy is implemented 
in the form of an element. In that work as is the case in this study, the interdomain 
interface boundary is discretized with a mesh of evenly-spaced pseudo-nodes (open 
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circles in Figure 2.6) that need not be coincident with any of the interface nodes (fdled 
circles in the figure) of any of the subdomains. 



The generalized element equations may be obtained by introducing the continuity 
requirements into the weighted residual statement. Eq. (2.36) can be rewritten over an 
element domain as 
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Note that in the potential energy formulation , the continuity of the secondary variables 


was satisfied through the subsidiary conditions obtained through the minimization of the 
potential energy. In this weighted residual formulation, the continuity of the secondary 
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variables is satisfied in a weighted residual sense and the Lagrange multipliers, and A , , 
are represented by weight functions in the form of the secondary and primary variables, 
respectively. 

The form of the equations for the finite element and finite difference applications 
differs by the form of the element shape functions and the approximation selected for the 
weight functions, <I> . For the generalized element expansion of subdomain i, the 
independent approximations for the element generalized primary variables, ( i.e ., 
displacements or velocities), interface secondary variables (i.e., fractions or fluxes), the 
weight functions associated with the secondary and primary variables, and the interface 
variables, are, respectively 

= N ; -u e . ; q, = R/ w ; ; Xt=R t ; i = T andv = T Ul 

where a is a vector of unknown coefficients associated with the secondary variables, q , 
and N, R, and T are matrices of interpolation functions for the element primary and 
secondary variables, and primary variables along the interface, respectively. The 
interpolation functions in the matrix R are assumed to be constants for linear elements 
and linear functions for quadratic elements. The interpolation functions in the matrix T 
are cubic spline functions. The derivation of this interpolation matrix is given in 
Appendix B, and the derivation of the geometry representation, T 1 , is given in Appendix 
C. Substituting these approximations into Eq. (2.45) yields an integral equation in terms 
of the weight function, <t>, which is given by 
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where integration over the common subdomain boundary, P, is considered only for 


element edges along that boundary. 
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Note that all of the element submatrices in the three-approximation formulation 
except for the k j matrix are identical to those obtained in the two-approximation 

formulation. The submatrix, k j , does not exist in the two-approximation formulation 

but is included in the three-field formulation. This submatrix is associated with the 
coupling of the primary variables along the subdomain interface boundaries to those 
along the interface. 

Assembling the element equations over the entire domain, enforcing continuity of 
the primary and secondary variables only within each subdomain and assembling the 
contributions along the element edges on the common subdomain boundary, and noting 
that u Cl and u e2 , and f ej and f e2 , andcq and a 2 are completely uncoupled, yields the 

system of equations given by 

Kj 
0 
0 

K Pi 

0 

or in a symbolic manner 

where K, u, and f are the assembled stiffness matrix, displacement vector and force 
vector for the entire structure, and K p , K s , Ki, m, and a are the assembled K p , K s , Kl, 
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ui, and a,- for all interfaces. 
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While it is convenient to represent the weighted residual form over the domain 
using a single equation, the system of equations, Eq. (2.46) is obtained from the 
individual weighted residual expressions over each of the subdomains and the constraint 
integrals. The first two matrix equations of the system of equations, Eq. (2.46) are 
derived from the weighted residual statement for subdomain i. That is, 


J k i 

A h 


d®; dUf d<t>i \ f t ^ 


- + - 


dx dx dy dy 


F>i - j ki dr 1 = lQ t O t da t + \q^i drf 


r 1 “i r / 

The third matrix equation of the system results from the reciprocity statement of the 
secondary variables. That is, 


J'fe+^ 2 )dr I =0 on r 1 . 
r 1 

The fourth and fifth matrix equations result from the continuity requirement for the 
primary variables, which is given by 

\X x {u Y -w^dr 1 =0 on r 1 
r 1 


jA 2 (^ 1 -i^dr 1 = 0 on r 1 

r 1 

For the finite element development, the weight functions are taken to be the finite 
element shape functions (i.e., <t> ; = N f ). For the finite difference development, the 
weight functions are taken to be the Dirac delta function (i.e., 

O; = (x - Xf , y - y t ) = (x,- ,>',•)). Thus, for i= 1,2, the finite element and finite 

difference stiffness matrices and force vector, k e , k p , k s , andf C; , for the three- 

approximation formulation are identical to those obtained for the two-approximation 
formulation for the respective discretization approaches. 
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Note that, for both of the discretization methods, the form of the coupling element 
matrices that are not in terms of the weight functions are independent of the method of 
discretization. That is, 

kp. =- jR/N.-dr 1 *, 

T e 

r 1 

and 

kj, = JT T R,- dT f 
r 1 * 

are of the same form for the finite element and finite difference discretizations. 

However, since the element shape functions, N,-, differ for the two methods, the interface 
matrices, k pj . , in general, are not identical. Moreover, in the finite element development, 

the weight functions, , are taken to be the finite element shape functions, N,-; thus, at 

T T 

the element level k s . =k , and at the global system level K s . =K 

i p i i P i 

The three-approximation derivation is more general as it allows for the coupling 
of the primary variables to an independent approximation. This attribute is particularly 
important in the heterogeneous discretization approach described in the next section. 

2.5.4. Multiple-Domain Modeling - Heterogeneous Discretization 

Heterogeneous discretization approaches make use of different discretization 
methods for at least two of the subdomains in which the domain is subdivided. There are 
many combinations of spatial modeling approaches; however, this work will focus on the 
coupling of the finite element and finite difference methods. 

Both the two- or three-approximation multifunctional formulations, discussed for 
the homogeneous discretization approach, are applicable to heterogeneous discretization. 
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However, as noted earlier, the three-approximation approach provides additional 
flexibility for the interface definition. Thus, only the three-approximation approach is 
presented. Hence, the multifunctional weighted residual formulation of Eq. (2.46) is 
used. Considering the two domains, upon which this discussion is based, one subdomain 
is discretized using the finite element method, and the other subdomain is discretized 
using the finite difference method. As before, for the finite element development, the 
weight functions are taken to be the finite element shape functions (i.e., <I>, = N ; -), and for 
the finite difference development, the weight functions are taken to be the Dirac delta 
function (i.e., <I>, = S j (x-x i ,y-yi) = S j (x i ,y i )). As expected, the set of element 

matrices becomes a hybrid of the matrices from the finite element method and the finite 
difference method. For completeness, these matrices are repeated here considering 
subdomain 1 as the finite element subdomain and subdomain 2 as the finite difference 
subdomain, and 


k e, = J h 

Q f 


3Nj r ^ 


dx dx dy dy 


|dDj and k e2 = & 2 


d 2 N T 2 


dx 2 
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k si =- W R 1 dpI and k s 2 =- R 2 * 


(2.47) 


and 


f e i = j N i T Gid^i + JlV^idr^ and f e2 = Q 2 (x i ,y i ) + q 2 (x h y i ), 

r? e 


and for the two domains, /= 1 ,2, 
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k p . =- JR/N.dr 1 * , 

r 1<2 

and 

ki. = Jt t r ; - dr 1 * . 

i e 

r 1 

2,6. COMPUTATIONAL IMPLICATIONS 

The two- and three-approximation multifunctional modeling approaches have 
been generalized such that they are applicable to both homogeneous and heterogeneous 
discretization approaches. Computational implications are presented in this section for 
the generalized system of equations, Eqs. (2.42) and (2.46). Implications specific to a 
discretization approach are highlighted, where appropriate. 

The assembled stiffness matrix K is a block diagonal matrix containing the 
stiffness matrices K, of each of the subdomains along its block diagonal. The interface 
“stiffness” matrix thus contains coupling terms that augment the stiffness matrices of the 
subdomains along the interface. The two- and three-approximation approaches yield 
systems of equations (see Eqs. (2.42) and (2.46)) of similar form and with the same 
attributes. Due to the use of Lagrange multipliers in the constraint conditions, the 
systems are neither banded nor positive definite. Therefore, standard Cholesky solvers 
can not be used, unless full pivoting is performed to obtain the solution. In addition, due 
to the generalization for the finite difference approximations, the system of equations is 
not necessarily symmetric due to different off-diagonal submatrices, K p and K s . The 
system unknowns in Eq. (2.46) consist of both primary and secondary variables given by 
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the potential function, u, and the secondary variable coefficients, a, respectively. 
Generally, the coupling matrices, K s , are of the order of the length of the interdomain 

boundary, which results in a marked difference in the magnitude of the off-diagonal 
terms of the system matrix compared to its diagonal terms. This characteristic produces 
an ill-conditioned matrix whose solution can cause difficulties for some general-purpose 
solvers. Hence, the coupling matrix should be scaled such that it is of the same order as 
the subdomain stiffness. The upper diagonal submatrix blocks contain uncoupled 
subdomain stiffness matrices. The symmetry of the subdomain matrix is determined by 
the choice of the weight function, O. For the finite element discretization, the subdomain 
matrices are symmetric. However, due to the elimination of fictitious nodes for the 
imposition of boundary conditions and loads in the finite difference discretization, the 
subdomain stiffness matrices, K,, generally are not symmetric, but they are positive 
definite and sparse. The coupling is accomplished through the introduction of the 
coupling terms in the matrices, K p/ and K s , for both approaches. The three- 

approximation approach requires an additional matrix, Kj . For the three-approximation 
approach, the number of additional degrees of freedom associated with the interface is 
generally small in comparison with the total number of degrees of freedom in the 
subdomains. Thus, modeling flexibility is provided at a relatively small computational 
expense. The computational expense in this study may be reduced additionally as the 
efficiency of new solution algorithms for the system of equations in Eqs. (2.42) and 
(2.46) is increased. 

The load transfer mechanism for finite element multiple-domain discretizations 

Of 

presented by Aminpour et al. is generalized for the multifunctional approach, herein. 
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This load transfer mechanism may be interrogated for the two- and three-approximation 
formulations by considering the first and second rows of Eqs. (2.42) and (2.46), 
respectively. For the three-approximation approach, the matrix equations of interest are 
given by 

Ki u i +K S] «i =fi 

K 2 u 2 +K S2 « 2 =f 2 

These equations can be partitioned such that they correspond only to the primary 
variables, u, on the interdomain boundaries. That is, u) represents a subset of u, ; 


hence. 


Kjui +K Sl «! =0 

K 2 u 2 +K S2 « 2 =0 


(2.48) 


where K, denotes interdomain boundary stiffness terms related to u, , and there are no 
forces on the interdomain boundary. The expressions given by the product term, K, u, , 
represent the internal fluxes at the 1 th interdomain boundary, and thus Eq. (2.48) may be 
written as 

?l = -K S] «j and f 2 =-K S2 « 2 . (2.49) 

For homogeneous discretization using the finite element method, substituting for 
K s ,. from Eq. (2.43) into Eq. (2.49) gives 


fj = - Jn/’Rj dT l6 «! =- jN/’q 1 dr I? 

r i e pb 

f 2 =- JnJr 2 dr I? « 2 =- jN 2 q 2 dr I<? 

T e je 


(2.50) 
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Examining Eq. (2.50) indicates that the evaluation of the internal fluxes is consistent with 
the evaluation of equivalent nodal fluxes in the presence of applied fluxes on the 
boundary. In addition, Eq. (2.50) substantiates that the secondary variables along the 
interface are represented by distributed fluxes for each of the subdomains. 

For homogeneous discretization using the finite difference method, substituting 
for K s ,. from Eq. (2.44) into Eq. (2.49) gives 

fl=-Rl«l=-qi (2.51) 

f 2 = -R 2 «2 = ~<l2 

Examining Eq. (2.51) indicates that the evaluation of the internal fluxes is consistent with 
nodal fluxes evaluated at points in the presence of applied fluxes on the boundary. In 
addition, Eq. (2.50) substantiates that the secondary variables along the interface for this 
approach are represented by nodal fluxes for each of the subdomains. 

For heterogeneous discretization using the combined finite element and finite 
difference methods, substituting for K s from Eq. (2.47) into Eq. (2.49) gives 

fj =- Jn/Ri dT 1 * «! = - jNfti! dr I<? (2.52) 

T e T e 

r 1 r 1 

f 2 = -R2«2 = - Q2 

Examining Eq. (2.52) indicates for subdomain 1 (the finite element subdomain), that the 
evaluation of the internal fluxes is consistent with the evaluation of equivalent nodal 
fluxes in the presence of applied fluxes on the boundary. Meanwhile, for subdomain 2 
(the finite difference subdomain), the evaluation of the internal fluxes is consistent with 
nodal fluxes evaluated at points. This reveals that for this multiple-domain approach, the 
secondary variables along the interface for subdomain 1 are represented by distributed 
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fluxes, while for subdomain 2 the secondary variables along the interface are represented 
by nodal fluxes. Thus, for this heterogeneous modeling approach, it is required to 
transform the interface secondary variables into equivalent quantities. 

2,7. VERIFICATION TEST CASE 

In this section, the multifunctional methodology for the scalar-field problem is 
demonstrated on a verification test case. The application is described, and the associated 
results and salient features are discussed. This application is considered a patch test for 
the formulation and verifies the applicability of the method for a configuration for which 
the solutions are known. Finite difference and finite element solutions for single- and 
multiple-domain configurations are presented to provide benchmark solutions for the 
multifunctional approach using homogeneous and heterogeneous discretization. 
Representative applications from the field of engineering science are presented in 
Chapter IV. 

2,7.1. Patch Test Problems 

The patch test has proven to be a useful discriminator of the convergence 
properties of finite elements and other discretization approaches. A patch test refers to 
any problem with an exact solution as a constant state for which the approximating 
primary variable is capable of reproducing. The fundamental concept of the patch test for 
the scalar-field problem herein is to subject a domain to boundary conditions that 
engender a linear or quadratic primary variable field and a constant or linear secondary 
variable field throughout the domain. For the governing differential equation of the form 
of Eq. (2.1), boundary conditions that serve this purpose are: 

i. Specified primary variable onf^ which emanate from a linear field as 
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u = a^x + a 2 y + 

or quadratic field as 
(l 2 ^ 

u = ai\x -y }+ a 2 x + a^y + a,Q 
where a v a v a,,. and a (j are arbitrary constants. 

ii. Constant or linear secondary variable on T' s 
q = b x x + b 2 y + b Q 

Given these boundary conditions, a solution is sought to the Laplace’s equation. This 
governing equation is applicable to a variety of problems in engineering science. For 
example, consider the solution for the primary variable, u(x,y), in a rectangular domain 
(see Figure 2.7) with boundary conditions of the forms indicated which yield the exact 
solution. 



Figure 2.7. Two-Dimensional Rectangular Domain. 

The problem is given by 

d 2 u d^u 
dx 2 + dy 2 


0 <x<a. 


0 <y <b 
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which is known as Laplace’s equation for a planar domain. 

Results of the analyses performed have been compared to appropriate reference 
solutions and are summarized in Table 2.2 using normalized values. A value of unity 
implies perfect agreement with the reference solution. Specified boundary conditions 
representing linear, bilinear, and quadratic potential functions are applied to the square 
domain. For all cases, the reference solution is the exact solution. For the linear case, a 
specified boundary condition of the form 

tt(0,y) = 2, u(a,y ) = a + 2, and q n (x,0) = q n (x,b) = 0 
has been imposed. For the bilinear case, a specified boundary condition of the form 
u(0 ,y) = y, u(a,y) = a + y, and q n {x,®)=-\ and, q n (x,b) = 1 

has been imposed. For the quadratic case, a specified boundary condition of the form 

2 2 2 

«(0,y) = -y , and u(a,y) = a —y , q n (x,0) = 0, and q n (x,b) = -2b 

has been imposed. Several analyses have been performed namely, (1) two single-domain 
analyses with individual finite element and finite difference discretizations, respectively, 
(2) two multiple-domain analyses with homogeneous modeling with individual finite 
element and finite difference discretizations, respectively, and (3) one multiple-domain 
analysis with heterogeneous modeling with combined finite element and finite difference 
discretizations. Results from these analyses are summarized in Table 2.2. In this work, a 
five-node central difference template and four-node quadrilateral finite elements are used 
to form the models. Spatial modeling is used consistent with single-domain modeling 
approaches with a (5 x 5) mesh and a (9 x 9) mesh. The coarse and fine models, shown 
in Figure 2.8, are used in the finite element homogeneous modeling. For the finite 
difference homogeneous modeling and the heterogeneous modeling, a finite difference 
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mesh is used that has the same number of mesh points as the finite element mesh in the 
respective domain. 



Figure 2.8. Spatial Discretization for Two-Dimensional Rectangular Domain. 

For boundary conditions consistent with linear and bilinear potential functions, 
the computed potential and flux results are exact for all analysis types. For boundary 
conditions consistent with a quadratic potential function, the error in the computed 
potential and flux is approximately 3% for the multiple-domain homogeneous finite 
element (MDFE) spatial modeling, and the error is approximately 1% for the multiple- 
domain heterogeneous modeling (MD/HM) with finite difference and finite element 
discretization. For the given boundary conditions and element configuration ( i.e ., square 
or rectangular elements), the single-domain finite element (SD/FE) model reproduces the 
exact solution using the bilinear finite element. However, for a general element 
orientation (i.e., quadrilateral elements), the bilinear element used does not reproduce the 
exact solution. Moreover, for the multiple-domain analysis, error is introduced when 
combining finite element models of different discretization along the boundary. This 
error is due to the use of a higher-order interpolation function (i.e., cubic spline) on the 
interface than that used to represent the potential on the finite element edges. The error 
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obtained using the heterogeneous model is smaller than that obtained for the 
homogeneous finite element model. This attribute is due to the ability of the finite 
difference model to represent accurately the potential function on the interface based on 
the higher-order shape function used in the generalization of the finite difference method. 


Table 2.2. Results of the Multifunctional Approach for the Patch Test Problems. 


Analysis 

Normalized Potential Function, u 

Normalized Flux, q x 

Order of Potential function 

Order of Potential 

"unction 

Type' 

Linear 

Bilinear 

Quadratic 

Linear 

Bilinear 

Quadratic 

SD/FE 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

SD/FD 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

MD/FE 

1.0 

1.0 

1.03125 | 

1.0 

1.0 

1.03125 

MD/FD 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

MD/HM 

1.0 

1.0 

.98958 

1.0 

1.0 

.98958 


SD/FE: Single-Domain with Finite Element discretization 

SD/FD: Single-Domain with Finite Difference discretization 

MD/FE: Multiple- Domain with Finite Element discretization 

MD/FD: Multiple- Domain with Finite Difference discretization 

MD/HM: Multiple-Domain with Heterogeneous Modeling (combined finite 

difference and finite element discretizations) 
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CHAPTER III 

MULTIFUNCTIONAL APPROACH FOR VECTOR-FIELD 

PROBLEMS 

3.1. GENERAL 

While a scalar-field problem is one in which the dependent variable is a scalar and 
requires only the specification of magnitude for a complete description, a vector-field 
problem is one in which the dependent variable is a vector of components and requires 
the specification of magnitude and direction. Many of the concepts outlined for the 
scalar-field problem in the previous chapter are readily extendable to the vector-field 
problem, which allows further generalization of the multifunctional approach developed 
herein. A representative example of the vector-field differential equation in two 
dimensions is considered, and the mathematical statement is formulated. The concepts 
developed here are directly applicable to one-, two-, and three-dimensional vector-field 
problems; however, only the two-dimensional development is included in the interest of 
brevity. The general form of the differential equation describing the vector-field problem 
governing the motion of a continuum is given by the equilibrium equation 

pb + VT = -^! (3.1) 

dt 

where the variables p, b, T and v are the material mass density, the body force per unit 
volume, the stress tensor and the velocity vector, respectively. Eq. (3.1) is subject to the 
natural boundary condition, t = T • n = t on U, and essential boundary conditions, 

u = u,onr p where the normal vector to the boundary T is given by n = n x i + n v j , and 
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n x and n y are direction cosines of the unit normals, i and j . In addition, t , and u are 
applied tractions, and prescribed displacements, respectively, and v is the initial velocity 
vector. The equilibrium equations must be satisfied within the domain. Note that instead 
of prescribing the tractions on the boundary, boundary conditions may be given in terms 
of displacement or velocity components. Furthermore, boundary conditions on r may be 
mixed (i.e., surface forces, t, maybe prescribed on one part of the boundary and 
displacements or velocities may be prescribed on another). The equilibrium equation and 
other governing equations of continuum mechanics are discussed in more detail in the 
following section. 

3.2, CONTINUUM MECHANICS FOUNDATIONS 

The conservation of mass, linear momentum, angular momentum, energy, and 
entropy give rise to field equations that govern the deformation and motion of a 
continuum, and these equations are given in the form of integral or differential equations. 
In deriving the governing equations, the starting point is a statement of the conservation 
principle applied to a “control volume” to develop the integral form of the equation and 
extract the differential form by using the divergence theorem. 

3.2.1. Principle of Conservation of Mass 

The principle of conservation of mass states that when the total mass of the body 
is unchanged for an arbitrarily small neighborhood of each material point, the mass is 
considered to be conserved locally. Hence, the rate of increase of the mass inside the 
control volume is equal to the net inflow of mass through the control surface. 
Mathematically, this principle is given by 
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A 


jff + v- W 


FV 


dV = 0 


Since the integral is equal to zero for arbitrary respective volumes, V, the integrand must 
be equal to identically zero everywhere in the domain. The resulting equation, known as 
the continuity equation, is well known in fluid dynamics and is given in a conservative 
form by 


+ V • (pv) = 0 or ^ + ^il = 0 
dt dt dxf 


(3.2) 


The differential equation takes on a slightly different form when the derivatives of 
products are expanded and the definition of the material derivative is considered. The 
resulting non-conservative form is given by 

dp „ A d/9 d\’j 

dt d t dX( 

If the material is incompressible so that the density in the neighborhood of each material 
particle remains constant as it moves, the continuity equation takes the simpler form 


V • v = 0 or — - = 0 (3.3) 

«*/ 

This is known as the condition of incompressibility, which is important in classical 
hydrodynamics and plasticity theories. The continuity equation is an important partial 
differential equation in all branches of continuum mechanics and the discipline-specific 
aspects are discussed in the next section. 

3.2.2. Conservation of Linear Momentum 

The equations of motion, valid in all branches of mechanics, are partial 
differential equations derived from the momentum principles of a collection of particles. 
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In this case, it is easier to use integrals over a given mass of material (the material 
volume, V') rather than over a given spatial volume (the control volume, V). The 
Reynolds transport theorem is used to replace the material volume with the control 
volume. The conservative form of this theorem is given by 


_d 
d t 


\ppdV 

V' 


— f<ppdV + j</>pvndS 
dt V S 


where 0 is the continuum property per unit mass and j p\ ■ n d.S' is recognized as the mass 

S 


flux. The conservation of linear momentum represents Newton’s second law and governs 
the motion of the continuum under the influence of the external effects. This principle 
states that the time rate of change of momentum is equal to the resultant force, F, acting 

on the body. Thus, F = where F is the resultant of all external forces and is given 

dr 

acting on a material volume as F = J pb dV' + |t dS ' , and L is the linear momentum 

V S' 

vector on a material volume given by L = J xp dV ' . First, expressing the conservation 

V 


of linear momentum over the material volume and then using the Reynolds transport 
theorem to express the equation in terms of the control volume yields the integral 
conservative momentum equation given by 

J pb dV + tft dS* = — J \p dV + <f \p\ ■ n dS* 

V S V S 


Using the divergence theorem and Cauchy’s formula, the conservative differential form 
may be obtained as 


pb + VT 


djpv) 


+ V • (pvv) 


at 



The non-conservative form of the differential equations is obtained by expanding the 
divergence operator, V • (pvv) , and making use of the continuity equation, Eq. (3.3), 
yielding 
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pb + V-T 


d(pv) 

d t 


or 


Pty + 



dx j 


d(pv ; ) 

dt 


3.2.3. Conservation of Angular Momentum 

The principle of conservation of angular momentum is used to show symmetry of 
the stress tensor, which is used to describe the state of stress of the continuum. In a 
collection of particles whose interactions are equal, opposite and collinear forces, the 
time rate of change of the total moment of momentum for the given collection of particles 
is equal to the vector sum of the moments of the external forces acting on the system. In 
the absence of distributed couples, the same principle for a continuum is postulated. 

Thus, 

J(rxt)dS* + J(rxpb)dF = — J(rx pv) dV 

S V 

where x denotes the vector cross-product operation. Upon expressing the cross products 
in indicial notation, transforming the surface integral to a volume integral (using the 
divergence theorem), and using the expression for the material derivative of a volume 
integral, the moment of momentum equation is reduced to 

^krs^sr ~ 0 

at each point where e^s is the permutation operator. This yields 


For z-1 


Tyi - T23 ~ 0 
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Forr=2 T^\ - 7n = 0 

Forr=3 Tn - Ti\ = 0 

establishing the symmetry of the stress tensor in general without any assumption of 
equilibrium or of uniformity of the stress distribution. However, the balance of the 
couple stresses is assumed. In reference 39, a proof is given for symmetry of the stress 
tensor involving the condition that the rates of change of the components of stress remain 
finite. 

3.2.4. Conservation of Energy 

The principle of conservation of energy states that energy is conserved if the time 
rate of change of the kinetic and internal energy is equal to the sum of the rate of work of 
the external forces and all the other energies that enter or leave the body per unit time. 
Such energies supplied may include thermal energy, chemical energy, or electromagnetic 
energy. Herein, only mechanical and thermal energies are considered, and the energy 
principle takes the form of the well-known first law of thermodynamics. Since the 
energy equation involves an additional unknown quantity, the internal energy, the 
equation is a useful addition to the equations of continuum mechanics only when it is 
possible to relate the internal energy to the other state variables; in traditional 
thermodynamics an equation of state furnishes the required relation. The first law of 
thermodynamics applied to a material volume may be written as 
K + U = W + Q 

where the superscripted dot, ('), represents the derivative with respect to time, and K is 

the rate of increase of the kinetic energy of the material volume , U is the rate of increase 
of the internal energy of the material volume, W is the rate of work done by the external 



forces on the material volume, and Q is the rate of heat added to the material volume. 


The individual variables are defined as follows: 


K = — f -pv • v dV' 
dt f 2 


U = — \pudV' 

dt yf 


W= jp bvdV'+ jTv-ndS" 
V' S’ 


Q = - jq • n dS 1 ' + J pr dV' 

S’ V 


where u is the specific internal energy, q is the heat flux vector and r is the radiative heat 
transfer per unit mass. Upon using Reynolds transport theorem to convert the material 
volume to the control volume and the divergence theorem to convert the surface integrals 
to volume integrals, and performing further algebraic manipulation, the energy takes the 
form 

Afi 

\p— dV = -jVqdV+ jprdV+ jT:DdR 

V dt V V V 


where the stress power, T :D , is the scalar product of the stress tensor, T, and the rate of 
deformation tensor, D. The differential forms are given by 



If only mechanical quantities are considered, the principle of conservation of 
energy for the continuum may be derived directly from the equation of motion. This 
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equation, referred to as the conservation of mechanical energy, states that the rate of 
increase of the internal energy equals the heat added per unit time plus the stress power 
that is not contributing to the kinetic energy. The equation is given by 

du 


p — = -V q + pr + T:D 
At 


(3.4) 


3.2.5. Second Law of Thermodynamics 

The second law of thermodynamics is automatically satisfied and includes the 
change in entropy of the continuum. The entropy is regarded as a measure of change of 
energy dissipation with respect to temperature. The relationship expressing conversion of 
heat and work into kinetic and internal energies during a thermodynamic process is set 
forth in the energy equation. The first law, however, leaves unanswered the question of 
the extent to which the conversion process is reversible or irreversible. The basic 
criterion for irreversibility is given by the second law of thermodynamics through the 
statement on the limitations of entropy production. For a general process, the energy 
equation and the second law of thermodynamics are combined yielding 

d? _ 1 dq Q 
At T At pT 

where dv is the change in the entropy per unit mass, T is the absolute temperature, — is 

dr 

the heat transferred per unit time per unit mass, 0 is the dissipative function obtained 
from 0 = T® Dy using the dissipative or deviatoric stress tensor T°, and the notation d 
is used to indicate that the quantity is not an exact differential. The deviatoric stress 
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tensor is defined by Ty = Ty - pSy where - p is the hydrostatic pressure. For a general 
process, Q > 0 

d? > 1 Aq 
At T At 

and for an adiabatic process, 



At 


where in each of the above equations, the equality condition holds for a reversible 
process and the inequality condition holds for an irreversible process. 

The general principles of continuum mechanics have been outlined in this section 
to provide a foundation for the basic equations governing the motion of general continua. 
In the derivation of the balance laws, no differentiation has been made between various 
types of substances. The character of the material is brought into the formulation through 
appropriate constitutive equations for each material with the constitutive variables being 
restricted in their regions of definitions. These and other discipline-specific attributes are 
outlined in the following section. 

3.3. DISCIPLINE SPECIFICS 

The constitutive equations characterize the individual material and its reactions to 
applied loads. Hence, in the following section, the discipline-specific attributes of solid 
and fluid continua and their impact on the general principles of continuum mechanics are 
reviewed. In addition, other salient characteristics of the governing equations for solids 
and fluids are discussed. 

All constitutive equations must be consistent with the general principles of 
continuum mechanics. While impact of the constitution of the continua is discussed for 
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all of the balance laws, emphasis is placed on the principle of conservation of linear 
momentum. This principle is the basis for the governing equations of the multifunctional 
approach presented herein. This law states that the sum of the body forces together with 
the sum of the contact forces is equal to the change of the linear momentum of the 
material. The law is used as the basis for describing the motion in both solid mechanics 
and fluid mechanics. 

3.3.1. Solid Mechanics 

The field of solid mechanics has traditionally been characterized by well- 
formulated analysis of mechanical phenomena occurring in engineering systems, 
combined with experiments that explore the basic concepts 40 . Herein, elasticity theory is 
the primary field of solid mechanics discussed. In classical linear elasticity theory, it is 
assumed that displacements and displacement gradients are sufficiently small such that 
no distinction need be made between the Lagrangian and Eulerian descriptions. It is 
further assumed that the deformation processes are adiabatic (no heat loss or gain) and 
isothermal (constant temperature). The conservation of mass states that the mass of a 
deformed piece of material is the same as the mass of the undeformed material. In 
elasticity, based on the small strain assumption, the density, p, in the deformed state may 
be approximated by the density, p 0 , in the undeformed state, and the conservation of mass 
is identically satisfied. 

Moreover, it is convenient to identify a material particle of the continuous body 
by giving its initial coordinates. The position coordinates, x, y, z appearing in the partial 
derivatives and the integrals in the foregoing derivatives are, however, the instantaneous 
positions. For an elastic body in equilibrium, they represent the coordinates of a particle 
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in its new position in the deformed body. When the strains and displacements are small, 
it may be possible that the equilibrium conditions are satisfied in the undeformed 
configuration of the body. The equilibrium differential equations are strictly applicable 
and the stress tensor is strictly symmetric for the nonpolar case only when defined in the 
instantaneous deformed position. Even in small strain theory of elasticity, it is necessary 
to take account of this attribute in applications where the instability may occur, as in the 
buckling of a column or a shell. Asymmetry of the stress tensor also occurs when there is 
distributed couple stress 6 . 

In ideal elasticity, heat transfer is considered insignificant, and all of the input 
work is assumed to be converted into internal energy in the form of recoverable stored 
elastic strain energy, which can be recovered as work when the body is unloaded. In 
general, however, the major part of the input work into a deforming material is not 
recoverable energy stored, but dissipated by the deformation process, causing an increase 
in the body’s temperature and eventually being conducted away as heat. When thermal 
effects are neglected, the energy balance equation may be written as 

dw 1 r n = ]_ T p 

dt p ij ij p ij ij 


The internal energy, u , in this case is purely mechanical and is called the strain energy 
density (per unit mass) 

d u = — Tiide a 

p J J 

A material body is said to be ideally elastic when the body recovers (under 
isothermal conditions) its original form completely upon removal of the forces causing 
deformation, and there is a one-to-one relationship between the state of stress and state of 
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strain. The generalized Hooke’s law relates the nine components of stress to the nine 
components of strain 

a \ ij = C ijkl e kl 

Symmetry of stress and strain reduces the number of material constants in the fourth- 
order tensor, C , from 8 1 to 36. The existence of the strain energy density functional 

further reduces the number of constants to 2 1 . The existence of three mutually 
orthogonal planes of symmetry reduces the number of constants to nine. Isotropy reduces 
the number of constants to two. 

For this special case, Hooke’s law reduces to 

cjy = [A SijS M + Ju(S ik Sj} + 8 n 8j k )]% (3.5) 

where 

_ r _ E vE 

2(1 +v) ’ (l+vXl-2v) 

For i=j= 1, the second and third terms of Eq. (3.5) are nonzero if k=l and 1=1. Thus, 

fjjj = Ae + 2ii£\i 

where e = £\ \ + £22 + £33 • For i = 1 and j= 2, the second term of Eq. (3.5) is nonzero if 
&=1 and 1=2 and the third term is nonzero if k= 2 and 1=1. Thus, 

<*12 = M £\ 2+^21 = 2 mi - 
Similarly, other components of stress may be defined. 

Noting that the linear strain-displacement relationship is given by 


e ij 2 ^ Ui,J + Uj,i ^ 



88 


One method of solution of the problems of elasticity is to eliminate the stress components 
in the equilibrium equations given in indicial notation as 

°ij,j+p{bi -V/) = 0, 

and using Hooke’s law to express the strain components in terms of the displacements. 
Eq. (3.5) maybe written, with no loss of generality, as 

® ij ~ /"O/ "h 

Solving the boundary-value problem involving 15 equations for 15 unknowns is a 
formidable task. There are several ways of formulating the problem in terms of fewer 
unknowns and fewer equations. The most straightforward method is to obtain the 
stresses in terms of displacement gradients, and then substitute into the equilibrium 
equations to obtain three second-order partial differential equations for the three 
displacement components. Therefore, in terms of displacements, 

®ij ~ u j ,i ) 


and 


v t = u 


i 


dt 2 


Substituting these expressions into the equilibrium equation yields 


or 


*UI,V S V Ujjj P^Pi ® 


Uj,ij PPi ® 


Noting that / is a dummy index in the term Uiju The equation may be written as 

^ j 9 ji ^j\ij ) Pfoi ~ ) — 0 
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This leads to the field equations of Navier 

(T + p)u + pipi — iif ) = 0 


or 


p\' 2 u i +(A + p)u jji + p{b t -Mf)-0 


(3.6) 


The conditions for the static equilibrium of an elastic body are described by an elliptic 
system of nine partial differential equations for the displacements and stresses. 


3.3.2. Fluid Mechanics 

Fluids whose constitution is described by linear constitutive relations are called 
Newtonian fluids. The subject of Newtonian fluids is generally referred to as fluid 
mechanics, which encompasses widely diverse topics including, but not limited to, 
motion of airplanes and missiles through the atmosphere, the flow of liquids and gases 
through ducts, and the transfer of heat and mass by fluid motion. The constitutive 
equations for these fluids are given by 

a ij=~ P8 ij +C W D kl 

where P is the thermodynamic pressure and Da are the components of the rate of 
deformation tensor 

Dfci = 2 ^*,/ + v l,k) 

For isotropic fluids, the last term in the constitutive equations may be written as 

C ijkl D kl = ^D rr Sy +2 pDy 


or 

6//A7 D/{l = /\ 8 ij 8 kl ^ 8 jl + Sjj Sjj- 


Therefore, 
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®ij ~ Pb ij + A \&ij P^ik^jl l^jk^Pkl 

and by evaluating the Kronecka delta parameters, 

°ij = ~ PS V + AD kk S ij + 2 fJDij (3.7) 


This is the Navier-Poisson law for a Newtonian fluid. 

As in linear elasticity, substituting the constitutive equation into the equation of 
motion yields 


- PjSy + fa>k,ki S ij + P^tjj + v j,ij )+ P( b i -*i)- 0 


expanding gives 


or 


_ P,i + U + p)v j,ji + P v i,jj + P( b i - Vi) -0 


P^i ~ ~P,i + U + PY’jJi + P v i,jj + P b i 


or in vector form 


p— = -VP +(A + p)V(V • v)+ pV 2 v + pb 
dt 


Using the Stokes condition, A = p , the equations reduce to the Navier-Stokes 
equations and are given by 




(3.8) 


or 

p — = -VP + — V(V • v) + pV 2 v + pb (3.9) 

dt 3 

In this form, the difference between the Navier equations of solid mechanics, Eq. (3.6) 
and the Navier-Stokes equations of fluid mechanics, Eq. (3.8) or (3.9), can be readily 
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considered. In Navier-Stokes equations, there is not only an additional pressure term but 
also the equations are nonlinear; this can be seen by examining the acceleration, 

dvj ()\'j 

V[ = = — — + Vj- jVt, and from the products of the density, p, and the acceleration, v , 


present in the equation. Additional nonlinearities are evident in the continuity equation 
given by Vfc* = 0 (V-v=0). In the linear theory of elasticity, this situation does not occur 


. d 2 u 
since V; ~ — — 

^ 2 


and p is taken as a constant. The Navier-Stokes equations together with 


the continuity equation form a complete set of four equations and four unknowns: the 
pressure, P, and the three velocity components, v,-. 

For steady and low-speed flow of an incompressible fluid (V-v=Da=0), for 
constant p and by making use of the divergence-free condition in Eq. (3.8) or (3.9), the 
governing equations take the form 

D kk =0 


- + P v i,jj + P f) i ~ 0 

However, these equations, often referred to as Stokes equations, may be written for two- 
dimensions in the most general form without using the divergence-free condition to 
simplify the equations. In so doing, the physical form of the natural boundary conditions 
is preserved. The form of these equations is given by 



(3.10) 
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Fluids often behave as though they are inviscid or frictionless. Therefore, it is 
useful to investigate the dynamics of an ideal fluid that is incompressible and has zero 
viscosity. For frictionless flow of an incompressible fluid, the equations, called Euler’s 
equations, maybe obtained from the general Navier-Stokes equations. Since in a 
frictionless flow, there can be no shear stress present and the normal stress is the negative 
of the thermodynamic pressure, the equations of motion are 

pv, = ~P,i + Ph 


or 

p^ = -VP + pb 

at 

For a general fluid, the character ( e.g ., elliptic, hyperbolic, or parabolic) of these 
equations of motion is determined by the sign of the discriminant. The Navier-Stokes 
system of equations, in general, is considered as mixed elliptic, parabolic and hyperbolic 
equations. The system of time-dependent Navier-Stokes equations is essentially 
parabolic in time and space, although the continuity equation has a hyperbolic structure. 
Therefore, they are considered a parabolic hyperbolic system. For the same reason, the 
steady-state form of the Navier-Stokes equations leads to elliptic-hyperbolic properties. 
In addition, the classification of the differential equation changes with the flow 
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characteristics (i.e., subsonic, supersonic, or transonic), which may create great 
difficulties in solution where part of the flow is supersonic and part of it is subsonic. 

3.4. SINGLE-DOMAIN FORMULATION 

As in the scalar-field problem, methodology for the vector-field problem is 
presented formulating the general method of weighted residuals for a single domain. 
Consider the equilibrium equation governing the motion, u, of a continuum 

dv 

pb + V-T = p^~ 

ot 

or in indicial notation 

(jjj' j + p(bj —Vj) = 0 in Q for z, j = 1,2,3 (3.11) 

in a domain, Q, bounded by T. In this work, the equilibrium equations of Eq. (3.1 1) 
describe the motion of a three-dimensional continuum. Hence, the indices, i and j range 
from the value of unity to three (i.e., i,j = 1,2,3). This range will apply throughout this 
development unless otherwise specified. In general, the boundary, T, can have mixed 

boundary conditions with the primary variables, u, prescribed on T 1 and the secondary 

S 

variable, the traction, t, prescribed on the remaining part of the boundary, T . In solid 
mechanics, the six stress components will be some general functions of the components 
of the generalized displacement 

u T =[« v w 0 X Oy 0 Z \ 

where u, v, and w are translational components and 6 X , Q y and 0 Z are rotational 
components. In fluid mechanics, the stress components will be functions of the velocity 
vector 

u T =[vi v 2 v 3 ]. 
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which has similar components to those of the displacement vector. Thus, Eq. (3.11) can 
be considered as a general equation of the form of A(u) = 0 . 

The method of weighted residuals is applied to the vector-field problem in this 
chapter in the same way as for the scalar-field problem of Chapter II. Hence, an 
approximate solution, u , is used in expressing <5tjj through the use of stress-strain and 
strain-displacement (or stress-rate of strain) relations, then the differential equation, Eq. 
(3.1 1), will no longer be satisfied, and this lack of equality is a measure of the departure 
of u from the exact solution. The lack of equality is called the residual, R, and is written 
as 

R i = °ij,j + P {b i ~ V/ ) * 0 for i,j = 1 ,2,3 . 

The residual is orthogonalized by a set of weight functions, <E>, and may be written as 

J*/®/ dD= j(<7 VJ +p(bi -v, •))<!>,• dO = 0 (3.12) 

£2 Q. 

n 

where the approximate solution is given by u = . As defined before, the 

m = 1 

functions, Y , are trial functions, and a m are arbitrary coefficients. The trial functions 

satisfy the homogeneous part of the essential boundary conditions, while V F Q satisfies the 

nonhomogeneous part. Using the general weighted residual form outlined in Chapter II, 

JOA(5)dD + j ®B(5)dr = 0 . 
q r 

where the residual in the satisfaction of the boundary conditions is orthogonalized by a 
secondary weight function, O . For the system at hand, a vector quantity is sought and 
the differential equation is a simultaneous system of equations. Here, 
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A(u) = (Jjj^ j + p(b/ - Vj ) = 0 , and the essential and natural boundary conditions are 
represented by 


By (u) = u - u = 0 or Uj -Uj =0 on 


and 


^2 (u) = t - 1 = 0 or tj —tf =0 on P s , 

respectively. Therefore, considering the approximate solution, u , we may write the 
general integral form of the differential equation governing the continuum motion as 

]<S> i {a« J +p{b i -v i ))d&+ dT^+ jO i2 { ti -F,)dT s = 0 (3.13) 


Q 




Note that the approximate solution may be selected to satisfy the essential and the natural 
boundary conditions and thus the boundary integral equations in Eq. (3.13) are identically 
zero. In this formulation, we will presume that the essential boundary conditions, i.e., 


u - u = 0 or — ti( =0 on 

are automatically satisfied by the choice of the function, u . Therefore, Eq. (3.13) is 
rewritten as 

jO^CTyj+pibi -Vi))dQ+ -/,)dr y =0 (3.14) 

where 0,2 = O,- . 

In the formulation herein, the order of differentiation on the stress term in the 
integral equation, Eq. (3.14), is reduced to obtain the weak formulation. Recognizing that 
the stress components are functions of the primary variable, u, which is approximated by 
u . For simplicity, the subsequent development is presented in terms of u. Application 
of the divergence theorem to Eq. (3.14) yields 
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-\a ij ^ i ,j6Sl+^a ij n j )^ i &T + \p(b i -v i )^ i d£l+ Jfo -F,)®,dr s =0. (3.15) 

q r o y s 

Note that the domain boundary is presumed to consist of boundaries on which the 

primary variable is specified and boundaries on which the secondary variable is specified, 

and r = + r s . Therefore, the boundary integral on r, may be expressed as 

)<!>,. dr= j(<T, y n,> f dr' , + 

r y p r s 

In the method of weighted residuals, the weight functions, O, satisfy the homogeneous 
boundary conditions for the primary variable, and thus, 0=0 on T . Therefore, the 
boundary integral on r is identically zero and Eq. (3.15) may be rewritten as 

- 1 \{ (7 ij n j) ( ^i dr' + \p(bi -Vj)0 ; - dQ + {((• -F/)0 ; dr v =0. 

Q Y s Q y s 

Since the weight functions, O and O , are arbitrary, they may be chosen, without loss of 
generality, such that, O = -O , and using the Cauchy formula, t t = (Tyti i , 

- dQ + jpibi -v / )® / dQ+ \t t ®,dT J =0 (3.16) 

Q Q y s 

or 

J <T if ® 1 9j dD = J p(v t - b { )0 ; dQ + J t { ® ,dT‘ s . 

q q r * 

The integral fonn of Eq. (3.16) is given for a general continuum. If the weight functions, 
O, , are selected to be virtual displacements or velocities, Su, , then Eq. (3.16) is given by 

- f ; Ail + \ p(bf — V/ )5ui dQ + \ti <Sw/dr 5 =0. (3.17) 

o q r * 


The term Stij ,j can be expanded to 
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) 


where e and or are symmetric and skew-symmetric tensors, respectively. These tensors 
are given by 


ij 0 i u i,j + u j,i ) anc * 0) ij 0 i u i,j u j,i )• 


In solid mechanics, these tensors represent the linear infinitesimal strain-displacement 
and linear infinitesimal rotation tensors, respectively. In fluid mechanics, the tensors 
represent the linear infinitesimal strain-rate of deformation and vorticity tensors, 

respectively. Noting that a., is a symmetric tensor and that the product of a symmetric 

y 

tensor and a skew-symmetric tensor is zero, Eq. (3.17) maybe rewritten as 


- J (Tjj8s r . dQ + J p(bi - Vj )Siij dQ + J/ ; - c»qdr s =0 (3.18) 

Q £2 r s 

Eq. (3. 1 8) represents the principle of virtual work where the first integral term represents 
the internal virtual work, the second and third terms represents the external virtual work 
due to body forces, inertial forces and surface tractions. 

In the virtual work development, the term virtual work is loosely used for fluid 
mechanics and has been included here to highlight the similarities between solid and fluid 
mechanics. Variational techniques for perfect fluids, non-Newtonian fluids and general 
Navier-Stokes equations are discussed in Finlayson 7 . In this work, concentration is given 
to the general weighted residual equations, Eq. (3.16), and these equations form the basis 
of finite element approximations, which will be presented briefly in a subsequent section. 

Thus far, the single domain formulation has been developed for the vector-field 
problem focussing on the momentum equation, which is applicable to general continua. 
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However, the motion of a fluid is governed by the conservation laws of mass, momenta, 
and energy. In general, these equations consist of a set of coupled nonlinear, partial 
differential equations in terms of the velocity components, temperature, and pressure. 
When the Reynolds number for the flow is very low, the nonlinear terms due to inertial 
effects can be neglected, resulting in a linear boundary value problem. Such a flow is 
called Stokes flow 41 (see Eq. (3.10)). When temperature effects are not important, the 
energy equations are uncoupled from the momentum (i.e., Navier-Stokes) equations. 

Thus, for isothermal flows, only the Navier-Stokes, Eq. (3.8), and continuity, Eq. 
(3.2), need to be solved. Hence, an additional equation expressing the continuity 
condition is included in the weighted residual formulation. In the interest of 
completeness, the formulation herein is described using a Newtonian fluid. The laws 
governing the flow of Newtonian fluids were reviewed in Section 3.3.2 in which the 
equations were specialized to viscous fluids that are subject to the assumption of 
incompressibility. Under these conditions, the weighted residual statement of the 
equation of continuity, Eq. (3.3), is given by 

J» // 6dD = 0 (3.19) 

Q 

where the residual in the continuity condition is orthogonalized by the weight function, 

„ dvi 

<E> , and u { . = V u = — — . Hence, for fluid mechanics, both Eqs. (3.16) and (3.19) are 

’ J dx j 

the weighted residual statements required to approximate the continuum motion. While 
for solid mechanics, since the continuity condition, Eq. (3.3) and likewise Eq. (3.19) are 
automatically satisfied, Eq. (3.16) is the only weighted residual statement required to 


approximate the continuum motion. 



99 


3.5. MULTIPLE-DOMAIN FORMULATION 

As in the case of the scalar-field problem of Chapter II, the domain of the problem 
is subdivided into smaller subdomains. Consider the equilibrium equation governing the 
motion, u, of a continuum 

Oyj + p(bi - V() = 0 in Q for i,j = 1,2,3 (3.20) 

in the entire domain, Q, bounded by T. For simplicity, the multiple-domain formulation 
is presented for only two subdomains, Qi and (see Figure 2.3) with a single interface 
boundary. Independent approximations and weight functions are assumed in each of the 
subdomains and continuity conditions are used to provide for a continuous solution 
across the domain. Thus, Eq. (3.20) is satisfied in each subdomain, independently, i.e., 

ajy y + P\(h\^ - V ■ ) = 0 in Q| and + Pi (bj 2 - - v^) = 0 in fi 2 

subject to the boundary conditions on the subdomain boundaries, r 1 and lA, and the 
superscripted numbers enclosed by parentheses denote the subdomain. In general, the 
boundaries can have mixed boundary conditions with the primary variable, u, prescribed 

on and the secondary variable, the traction, t, prescribed on r . These boundary 
conditions may be written as 

U| -U| = 0 or iSp =0 on T? and t — t j =0 or -t^ =0 on rj s 
and 

u 2 -U 2 = 0 or u =0 on T? and t - 1 2 = 0 or - 1 j 2 ^ =0 on fj. 
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For the multiple domain case, the boundary at the interface between the two subdomains 
is denoted T 1 . Hence, the subdomain boundaries, T*, are presumed to include three 
boundary types, and these boundaries are given by 

n =r t +r t +r t for k = 1,2. 

Here, the boundary on the subdomain common boundary is assumed to represent the 
same geometry and thus, = T 1 . The residual for each domain is ortho gonalized by a 

set of weight functions, and is written as 

Q| 

and 

j(aM + A(*P ) -ff ) >M 2)d n 2 =° 

S! 2 


~ n ~ n 

where the approximate solution is given by ii\ = Y. a \m^\m an d #2 = Y^ a 2rn^2m ■ The 

m = 1 m = 1 

functions, T lm and 4' 2m , are the trial functions, and a {m and a 2m are sets of arbitrary 

coefficients. Using the general form outlined previously, ( i.e 

JOA(S)dfl + fOB( u)dr = 0 ), for each subdomain, one may write 
Q r 

J<J>^A(#))dQ + =0 for k = 1,2 

^ k ^ 

Therefore, considering the approximate solutions, u^and u' 2 ^ , the general integral form 
of the differential equation governing the motion for subdomain 1 is given by 
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J ®! 1 ) ( 4 ) y + A (#-#))* 2 1 + J 

+ J : ')dTf = 0 


and for subdomain 2 as 


Q- 


+ 


1 W) 


r 2 


dr 2 =0 


r 2 


Again, the essential boundary conditions, i.e., 

rW _»(!) = i 


and 


uj - uj = 0 or w ; - -Mj- =0 on Tf 


ii 2 - u ; = 0 or uj 2 ^ -uj- 2 ^ = 0 on 


are identically satisfied by the choice of the functions, uj and u 2 . Therefore, for 
subdomain 1, Eq. (3.21) is rewritten as 

+A(*/ (l) -V/ (l) ))dOi + J vPifP -iP)ar f =0 




r f 


where ojV = of- . Similarly, for subdomain 2, 

j®i 2) (4H + P 2(*, (2) -> > ( 2) ))< i£i 2 + j =° 


^2 


1 2 


where 


(3.21) 


(3.22) 


(3.23) 


(3.24) 
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The order of differentiation on the primary variable in the integral equations, Eq. 
(3.23) and (3.24), is reduced to obtain the weak formulation. Using the divergence 
theorem Eq. (3.23) can be rewritten, for subdomain 1, as 


- J dDj + <f(c7^n^)o^dr 1 + -v^)<D^ dDj 

rj £2j 


(3.25) 


and similarly, for subdomain 2, 


- J o[ 2 ]cr ; 5 2) dQ 2 + ^fnf^ dr 2 + \p 1 {b\ 2) -v^®f ) dD : 
’ r 2 


+ 


J op’fT’-'f’Ki =' 


(3.26) 


Recall that the boundary T is presumed to consist of boundaries on which the primary 
variable is specified and of boundaries on which the secondary variable is specified, and 


boundaries at the subdomain interface, and for subdomain k, IT = T? + rf + f 


I 


Therefore, the boundary integral on I* may be expressed as 



ri 


Noting that, <E>^. = 0 on YP . Therefore, the boundary integral on Tp is identically zero, 
and Eq. (3.25) can be rewritten, for subdomain 1, as 
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- j ®i‘M 1)da i + 1 (4‘M ,) M 1) <n?+ 1 (4M 1) M 1) 

Qi r * 5 r 1 


dr 1 


+ 




Oi 


Since the weight functions, and <t> > 1 ' , are arbitrary, they may chosen such that 


dr! 1 ) = -cftW 


<D> ' = -<f>) ’ , and using the Cauchy formula, t> ' = erf ' , 




- j ®W >d °i + 1 (4M 1) H 1)drI+ j 

«i r 1 Tj 

+ jAC^-v^oj^dDj =0 


(3.27) 


Qi 


Similarly, for subdomain, Q.%, 


- I <f u 4 2> dfl 2 + 
0,2 



dr 1 + j (p*®* 2 * cir| 


+ JP 2 (*i 2 ) -' > S 2 ) )®i 2) <1^2=0 
^2 


(3.28) 


In the two-approximation formulation for the scalar-field problem, the two 
primary field variables, ui and U2 are approximated independently, and continuity 
requirements between these two fields are satisfied at the subdomain interface boundary. 
The three-approximation approach, which makes use of a third approximation field for 
the primary variables along the subdomain interface boundary in addition to the 
approximations given along the boundary of the subdomains, is most general. Hence, 
only the three- approximation approach will be discussed for the vector-field problem. 
This primary variable, v, along the interface is assumed to be independent of the primary 
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variables, ui and 112 , of the subdomains to which it is attached. These independent 
approximations give rise to continuity requirements along the interface of the form 

v - uj = 0 or v z - - uf' = 0 on r 1 

V-U 2 = 0 or Vj -Wj- =0 on r 1 
These constraints can be satisfied in the integral sense as 

J ^(v - uj) dr 1 = 0 or { - “PW 1 = 0 on pI (3.29) 

n n 

Jx. 2 (v-u 2 )dr 1 =0 or j4 2 T< -M P) drI =0 on pI (3.30) 

n r i 

where A^and are Lagrange multipliers or weight functions in the form of the 

secondary variable along the interface. An additional continuity requirement in terms of 
the secondary variable along the common subdomain boundary is required. These 

secondary variables, and tj- 2 \ are assumed to be independent of each other. These 

independent approximations give rise to continuity requirements along the interface of 
the form 

f/ 1 )+f/ 2 )=0 on r 1 

These constraints can be satisfied in the integral sense as 

=0 on r 1 (3.31) 

r 1 

where is a Lagrange multiplier or weight function of the form of the primary variable 
along the interface. Combining Eqs. (3.27) and (3.28) for the entire domain, including 
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the three continuity integrals at the interdomain boundary, Eqs. (3.29), (3.30), and (3.31), 


and recognizing that and = [a^n'p 


yields 


- /(®S4° K - iW5H 2) - i'/'M" drI - j<, ( 2 M 2) dr ‘ 


Qj 

+ 


Q2 


j4 1) ( v i""i ,) )dr I + j4 2) ( v /-“( 2, ) dr ' + /ti(</ 1)+ 'i (2) )dl' 1 




+ 


(3.32) 


r l 


Jft(*/ (2, -«! 2) )®, (2) <W2 + << / (2) 4>S 2, dT| 

“2 II 

In addition, for fluid mechanics, the continuity equation is given and satisfied 
independently over each domain as 

u\ [ \ =0 in fii and u\ 2 ) = 0 in Q 2 
l J 1 hj Z 


The weighted residual statements over the domains are given by 

JwfJ.O^dDj =0 and |wj 2 )4>j 2 ^ dD 2 = 0 . (3.33) 

£2i Q 2 

Here, note that no integration by parts is used on the continuity equations, and no 
relaxation of the differentiability on u can be accomplished since the resulting boundary 
conditions would not be physical. Combining Eqs. (3.33) yields 

dflj + dd 2 = 0 . (3.34) 

£2j Q 2 

The integral form of Eq. (3.32) forms the basis of finite element approximations for solid 
mechanics, and both Eqs. (3.32) and (3.34) form the basis for fluid mechanics. These 
finite element approximations as well as other approximations will be discussed in more 


detail in the next section. 
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3.6. SPATIAL MODELING FOR MULTIPLE DOMAINS 

Spatial modeling for multiple domains using the finite element and finite 
difference methods for the approximation of the vector-field problem is outlined in this 
section. A brief overview of discretization methods is given followed by spatial 
modeling for solid and fluid mechanics domains. 

3.6.1. Overview of Discretization Methods 

Finite element and finite difference discretization methods for the vector-field 
problem are outlined in this subsection. For a more detailed discussion the reader should 
consult the literature. 

The finite element method 

The finite element method for the vector-field problem is developed in the same 
manner as for the scalar-field problem. In the vector-field problem, the dependent 
variable in the integral equations is a vector of components. In general, the inplane 
vector components ( e.g displacements parallel to the x and y axes) are approximated by 
the same shape functions. For isoparametric elements, this approximation is the same as 
that taken for the shape. For the elasticity problem, the consideration for the strain- 
displacement relation, the Jacobian transformation, and the displacement gradient 
interpolation results in a more complex (the product of three matrices) set of equations 
than for the scalar field. 

The finite difference method 

The finite difference method is ideal for solving the governing partial differential 
equations of a continuum. It represents a variety of equations in engineering science; 
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however, the method has not been used in solid mechanics to the same degree as the 
finite element method . The decline in the use of the finite difference method in solid 
mechanics is largely due to the limited flexibility of its treatment of boundary conditions. 
Most finite difference developments avoid the general problem of boundary conditions in 
one of the following ways: (1) a scalar problem, such as those of the previous chapter, is 
solved as an example and the boundary conditions are incorporated in the analysis using 
arguments based on symmetry of the independent variables in the derivative 
approximations or (2) an example is chosen with fixed boundaries to eliminate the 
presence of fictitious points. The lack of an intuitive procedure for elimination of the 
fictitious or external grid points introduced when a central difference operator is applied 
to a boundary point is one cause of the deficiency in the method. For the vector-field 
problem discussed herein, a 3x3 central difference template is used to evaluate the 
momentum equation, Eq. (3.20). An approach for eliminating the fictitious points based 
on physical arguments is presented in reference 43. The fictitious nodes are replaced by 
boundary tractions using a set of constitutive equations and the primary variables in the 
continuum. These points can then be eliminated, and the boundary tractions are 
introduced into the finite difference model. An alternative approach is to construct 
special forms of the difference equations for grid points at or near the boundaries 44 . 

These forms make use of forward or backward difference operators to express differential 
forms. In general, standard forward or backward difference operators have higher-order 
truncation error than the central difference operators used for the differential equation. 
Hence, special forms using additional interior grid points are constructed such that the 
operators have the same order of truncation error as those operators used for the 
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differential equation. The latter approach is used in this work and will be discussed in 
some detail in the discussions of the patch test application given in this chapter. 

3.6.2, Overview of Single-Domain Spatial Modeling 

For a single domain, the finite element equations may be obtained by rewriting 
and manipulating slightly Eq. (3.16) over an element domain as 

JO /v -<7ydO e - J <I> iP(bj - Vj ) d£2 e - Jo^-dr^ =0 (3.35) 

Q e r s e 

where a., are the approximate stress fields produced by the stress-strain and strain- 
displacement (or rate of deformation) relations and approximating the primary variable 
over the element domain by u = Nu e . 

General finite element development 

Using the Galerkin method, the weight function is given by <l> = N . Substituting these 
approximations into the integral equation given in Eq. (3.35) and writing in matrix form 
yields 

j3N T <rdD e - jN T p(b-v)dD e - jN T tdr* e =0 (3.36) 

Q e r s e 


where 3 is the operator matrix defined, in general, by 
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d_ 

dx 


d = 


0 

0 

d 


dy 

0 

d_ 

_dz 


0 

d_ 

dy 

0 

_a_ 

dx 

d 

dz 

0 


0 

0 

_a_ 

dz 

0 

d_ 

dy 

d_ 

dx_ 


and the stress vector a is given by 

o = [<7n a 22 G 33 a \2 °23 a \3 F 

General finite difference development 

Recall that in the finite difference methods, derivatives are approximated by 
difference expressions that transform the derivatives and consequently the partial 
differential equation to algebraic expressions and equations, respectively. Upon 
substitution of the approximation function into the differential equation, the equations can 
be recast in weighted residual form by selecting <I>, = S(x - x { ,y - y { ) . Note that the 
subscript i on the weight function is used to denote the subdomain, while the subscript i 
on the coordinate values, x and y, is used to represent the point in the physical domain at 
which the Dirac delta function is evaluated. This nomenclature is used throughout the 
mathematical formulation presented here. The weighted form of the residual reduces to 
the evaluation of the partial differential equations using the approximate solution 
evaluated at the N selected mesh points 

For a single domain, as in the finite element method, the finite difference 
equations may be obtained by interrogating the weighted residual equations over an 
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element domain where the element, e, surrounds grid point i (see Figure 2.5). The 
approximate solution for the primary variable is given by 

M 

u = 'Yj^m u m oru = Nu e 
m = 1 

where Mis the number of shape functions over the element. The weight function, O, is 
given by the Dirac delta function, S\x - x t ,y - y t ) = f>(x ; -,_y ; - ) . Therefore, Eq. (3.35) 
becomes 

J 8{xt,yi \ j <*ij ~ J 8{xi,yi )p(b i - (\if 

n e Q e 

- \S(x i ,y i )t l dP s ' e =0 

r s e 

and upon making use of properties of the Dirac delta function, 

0g,j \x=xi - p[t>i Ui ,yt)- Vi (xj ,Vj )] - 1( (. x t ,y t ) = 0 . (3.37) 

y=y i 

This equation and the equations related to the finite difference formulation that follow are 
evaluated at point (X/,y ; ) where i denotes a point in the physical domain, and no 
summation is implied over the x ( - terms. Eqs. (3.36) and (3.37) are applicable to a general 
continuum irrespective of its physical constitution. Discipline-specific constitutive 
relations are considered at this point to continue with the finite element and finite 
difference developments specific to solid and fluid mechanics. Each of these 
developments will be discussed in turn. 

Solid mechanics - finite element discretization 

For solid mechanics, the constitutive relation relating stress and strain is given by 

c = E(c-c 0 ) + c 0 
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where the strain vector 


£ -[ e ll e 22 e 33 2e l2 2e 23 2 ^1 3 ] T » 

E is a matrix of material stiffnesses, and o fJ and £ 0 are initial stress and strain quantities, 

respectively. The strain-displacement relation is given by 

£ = 3u = 9Nu e = Bu e . 

Implicit in the definition of B is the use of the Jacobian matrix to transform from 
Cartesian coordinates to element natural coordinates used in the shape function 
development. In addition, in solid mechanics, the acceleration of the continuum is given 


by v = ii = 



Moreover, the second time derivative of the primary variable over the 


element domain is approximated by ii = Nii e . Substituting the stress-strain, strain- 
displacement relations and the acceleration into Eq. (3.36) yields 


r \ r 

Ti 


JB 1 EB 6£l e h e + 




Q e 


J/>N T N<K2* ti e = jB T E£ 0 dn e - jB T a 0 dQ £ 


Q. e 


Q e 


(3.38) 


+ jN T /*dQ e + jN T tdr* 


Q £ 


or 

k e u e +m e ii e =f e 

where k e is the element stiffness matrix, m, is the element mass matrix, 11, is the vector 
containing the generalized primary variables, ii e is the vector containing the second time 

derivative of the generalized primary variables, and f e is the element force vector 
containing the generalized secondary variables. Note that the acceleration term can be 
considered as an inertial force and included as part of the element force vector. 
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Assembling these element equations over the entire domain and enforcing continuity of 
the primary variable at the interelement boundaries yields the system of equations given 
by 

Mii + Ku = F 


YldctYl Tldcffl 

where K= £ jB T EBdQ e ;M= £ JpN T NdD e ; u is the assembly of all of the 
1 Cl e 1 Q e 


nodal degrees of freedom associated with the primary variables; ii is the assembly of all 
of the nodal degrees of freedom associated with time derivative of the primary variables. 


modes 


and F = X jB A E£ 0 dQ e - jB i c 0 dQ <? + JN A p(b -ii) dQ e + jN T tdP 


a e 




cr 


Solid mechanics - finite difference discretization 

For solid mechanics, making use of the stress-strain and strain-displacement 
relations, and substitution of the primary variable approximations into Eq. (3.37), the 
element equation becomes 


3 T EB 

x=x r - 

u e + 

1 

X 

II 

1 

T i T i 

n e = d E£q \x=Xi — 5 (Tq . 


y=yi _ 


L y=yi J 

y=yi 


+/Mxt>yt)+t(xt,yt) 


For the second derivative difference approximation, the number of shape functions, M= 3 
and u J ={«/_! u t u i+l }. 

Therefore, as in the finite element method the difference equations may be written in the 
form 


m e u e +k e u e =f e 

where k, and m e are the finite difference “element mass and stiffness” matrices, u e is the 
vector of generalized primary variables, ti e is the vector of time derivatives of the 
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generalized primary variables, and f e is the finite difference generalized force vector. 
Assembling the element equations yields 

Mii + Ku = F 


where 


* 

II 

Ml* 

s 

3 t eb 

X=Xf 

Nelem 

,M= x 

1 

_Z_ 

X 

II 

j 

, u and ii e contain all of the nodal 

1 


y=yt _ 

i 

X, 

II 

Ps 



degrees of freedom associated with the primary variables and its time derivative, and 
F = X d Es 0 |x=x/ -d c 0 1 X=X ( + /9b(x,- ,y t ) + t (xj , y, ) . 

i L y = >’i y = yt 

Fluid mechanics- finite element discretization 

For fluid mechanics, the constitutive relation relating stress and the rate of 
deformation, Eq. (3.7), for an incompressible fluid is given by 

where the viscous stress vector, r, is given by 

T = [tH t"22 7 33 r 12 r 23 t 13 ]> 

u denotes the velocity vector, P is the pressure, and I is the identity matrix. The viscous 
stress is given by t = 2/yD where is the shear viscosity of the fluid and D is the rate of 
deformation tensor whose components are given by 

D ij=\( u i,j +u j,i) ( 339 > 

Hence, the rate of deformation is related to the deformation and may be expressed in the 
same form as the strain-displacement relation as 
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D = 9fU = 9fNu e =BfU e 

where 9f is a differential operator defined by 9f = T9 and Bf = T3N . The 
transformation matrix, T, is used to introduce the scalar multiple of the shear components 
of the rate of deformation (see Eq. (3.39)) and is symbolically defined as 


T = 


1 0 0 0 0 0 

0 1 0 0 0 0 

0 0 1 0 0 0 

0 0 0 ^00 

2 


0 0 0 0 


I 0 


0 0 0 0 0 ^- 


In addition, in fluid mechanics, the acceleration of the continuum is given by 
dv 3v 

v = — = 1- v • Vv . Moreover, the time derivative of the primary variable over the 

dt dt 


element domain is approximated by v = u = Nu e and P = NP e . Substituting the 


constitutive and rate of deformation relations along with the acceleration into Eq. (3.36) 
and rearranging yields 

\ f \ 


( \ ( 
T tvt j r\S 


J pH 1 N dD e 
D. e 


M e + 


J /?N T (Nu e )N 6£l e 


Q e 


M e + 


j2//B 1 B f dD £> |u 

U. e 


f \ 

JB t NI d Q e 

Q. e 


(3.40) 


p e = jN T pbdQ e + jN T tdP s 


Q e 


or 

nigUg+CgUg+kgUg-qgPg =f e 

where the element matrices lq. and m e and the element force vector, f e , are of similar form 
as those obtained in the solid mechanics development, c e , is a nonlinear element matrix 
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resulting from the total derivative of the velocity, u e is the vector containing the 
generalized primary variables, P e is the vector containing the element pressure variables, 
and u e is the vector containing the time derivative of the generalized primary variables. 

Hence, in the fluid mechanics development, the rate of change of the velocity - u is 
analogous to the second derivative with respect to time of the displacement (ii in the 
solid mechanics development). Moreover, the first integral term of Eq. (3.40) can be 
thought of as an inertial force. Assembling these element equations over the entire 
domain and enforcing continuity of the primary variable at the interelement boundaries 
yields the system of equations given by 

Mii + C(u)u + Ku-QP=F (3.41) 

nelem _ nelem _ 

where K= £ jlju B 1 B f dT2 e ;M= £ { pN 1 NdD e ; 

1 Q. e 1 Ll e 

nelem nelem 

C= X JyON T (Nu e )NdQ e ; Q = £ JB T NI 6D. e ; u is the assembled vector 
1 Q e 1 Q e 

of all nodal degrees of freedom associated with the primary variables; P is the assembled 
vector of all nodal degrees of freedom associated with the pressure, u is the assembled 
vector of all nodal degrees of freedom associated with the time derivative of the primary 

modes T e 

variables and F = X {NVbdQ e + ^ l tdT s . 

1 Q e T s e 

In addition to the element equations for momentum, Eq. (3.40), the element 
equations for continuity must also be developed from Eq. (3.19). Using the Galerkin 

/V /V 

method, the weight function corresponding to the continuity equation is given by O = N . 
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Substituting the approximation for the weight function and the primary variable into Eq. 
(3.19) and writing the equation over the element yields 

( \ 


[B T Ndf2 e 


= 0 


Sf 


or 


qju e =0 


Assembling these element equations yields 

-Q T u = 0 (3.42) 

Equations (3.41) and (3.42) can be combined into one system of equations and written in 
matrix form as 


M 

0 

jui r 




0 

0 

W L 


C(u) + K 
/~»T 




or in a more symbolic form as 


MU + KU = F 


(3.43) 


where U = {uj U 2 U 3 P} T . Hence, the equations for fluid mechanics may be 
expressed in the same form as the equations for solid mechanics. Note that the system of 
equations, Eq. (3.43), is referred to as the primitive-variable model, the pressure- velocity 
model, or the mixed model . This mixed model results in a system that is nonpositive 
definite because of the zeros appearing on the main diagonal. In addition, the 
interpolation used for the pressure should be one order less than those that appear for the 
velocity field 41 . Furthermore, the pressure approximation may be discontinuous across 
interelement boundaries. In addition, because different orders of approximation are 
typically used for the velocity and pressure fields, the pressure may not appear at every 
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node of an element, which can complicate the assembly process. An alternative 
formulation, called the penalty function formulation’ 2,41 , circumvents this situation by 
treating the continuity equation as a constraint among the velocity components. This 
formulation is developed here for finite element discretization. 

From the weak form in Eq. (3.16), a functional describing the continuum motion 
can be obtained. The linear and bilinear forms of the functional over an element when 
the two-dimensional velocity field, (v h v 2 ), satisfies the continuity constraint, Eq. (3.19), 
is given by 


l( 0 1 , 0 2 )= \p(bi -vj)®! dQ e + jp(b 2 -v 2 )® 2 dO e + JfjOjdT^ 

n e Q. e r s e 

+ \t 2 o 2 dr s 


5((0i,0 2 ),(vi,v 2 )) = // J 2 

Q. e 

+ JU J 

Q e 


dOj dvj d0 2 dv 2 


dxj dxj dx 2 dx 2 


dQ £ 


"dOj 

dO 2 ' 
+ 


di’ 2 ' 
+ — — 

dxo 

9xi 

> 

[ dx 2 

dxj J 


dQ 6 


Note that the pressure does not appear explicitly in the bilinear form. The quadratic 
functional is given by 


7 ( v l >v 2 ) = , v 2 ), (vj , v 2 )) - L{v x , v 2 ) 


= f* I 

Q e 


dvj 

dxj 


+ 


dv x 

dxo 


1 

+ — 
2 


dv x _ + dv 2 _ 
dx 2 dxj 


diT 


- jp(b x -vj)vj df2 e - Jp(b 2 -v 2 )v 2 dQ e - Jt 1 v 1 dT^ 
Q e Q. e T s e 

~ \h V2dr s ' C 


(3.44) 



The equations governing the flow of viscous incompressible fluids, Eqs. (3.16), and 
(3.19), are equivalent to minimizing of Eq. (3.44) subject to the constraint 
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\ dvj dv 2 A 
G{v u v 2 ) = — + — = 0. 

0X2 (jX\ 

In the penalty function method, the constrained problem is reformulated as an 
unconstrained problem by minimizing the modified functional 

\2 


Im ( V 1 ’ v l)= /(vi , V 2 ) ■ + \ Ye | RV] , V 2 Y 

Q. e 


dD 


where the penalty parameter, y e , can be chosen for each element. The necessary 
conditions for the minimum of I m is 81 m = 0 or <?> V] I m = 0 and 8 Vl I m = 0 . 

where 8vi and §V 2 denote the first variation with respect to the velocity components, Vi 
and V 2 , respectively. Therefore, 


\ l m = I + Ve J G ’( v l > v 2 K, G(v { ,v 2 )df T 

Q. e 


= 1 


dSv 1 dvi 


d.X'j dxj dx 2 


ddvq f dx\ dv 2 x 
dx 2 r)xj 


d Q e 


J p(bi - 1 >! )<&! d Q e - J t x dv x dr- 5 


Q. e 

+7e I 


dSv/ 


Q 


e ^1 


dvi dv? 
— -+ — - 
dxq dx 2 


dQ = 0 


(3.45) 


and 
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8 vy I m = 8 v, J + Ye 1 G(v\ , v 2 ]S V . G(v\ , v 2 )dQ 


-v 2 

= J 

£2 e L 


O e 


3<5v 2 dv 2 3<5v 2 

2 // — - + JU 1 


3x 2 d x 2 9xj 


9v i 3v 2 

1 ■ + ■ 


3x 2 3x| 


dfr 


- Jy0(^2 -^ 2)^2 d^ e - J?2^2 drV 


(3.46) 


Q e 

+7e J 


98v 2 r 


Q 


e d “^2 


3vi 3 vo 

1 + ^ 


A 


3x| 3x 2 


dfl = 0 


These two statements, Eqs. (3.45) and (3.46), provide the weak forms for the penalty 
finite element model. While, the pressure does not appear in the weak forms explicitly, it 
is part of the boundary tractions, t\ and t 2 . The penalty finite element model is obtained 
using Eqs. (3.45) and (3.46), the approximations for the primary variable and the time 
rate of change of the primary variable, v = u = Nu e and v = u = Nu e , respectively, and 

by choosing dv\ = <5v 2 = N . Assembling these element equations over the entire domain 
and enforcing continuity of the primary variable at the interelement boundaries yields the 
system of equations given by 

Mu + C(u)u + Ku + Su = F (3.47) 

nelem nelem T 

where K= £ |2//B 1 B f dD e ; M= £ J/oPTNdQ*; 


1 


Q. e 


1 


Q e 


riplpm tipIpyh 

C= £ J pN T (Nu e )NdD e ; S = £ y\ N, X] T N, X2 d Q e ; u is the assembled 
1 Q e 1 Q e 

vector of all nodal degrees of freedom associated with the primary variables; u is the 
assembled vector of all nodal degrees of freedom associated with the time derivative of 
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modes _ T _ e 

the primary variables, F= £ jNpb 6£l e + JN 1 1 dr 15, , and N, x . denotes 

1 Q e r s e 

differentiation with respect to the independent variable, x t , i= 1, 2. 

Eq. (3.47) may be represented in a more symbolic form as 

MU + KU = F 

where K = C(u)u + K + S . Note that this penalty finite element method yields a system 
of equations in terms of the primary variables, u, and does not include the pressures, P. 
The pressures may be obtained from the computed velocity field by 



where ( V\y,V 2 y ) is the finite element solution of Eq. (3.47). 

Fluid mechanics - finite difference discretization 

For fluid mechanics, making use of the stress-rate of deformation constitutive 
relation, and substitution of the primary variable approximations along with the 
acceleration into Eq. (3.37), the element equation becomes 

/}N|x=x,- u e + yo(Nu e )N|x=x,- u e + 2//9 T TdN x=x . u e - 9 T IN X=X . P e 

L y=y t J L y=yi] [ y= ;.\ [ y i;.\ 

= Mxi>yi) + i(xi,yi) 

Also, considering continuity, 

9 T I N x=x . Ug = 0 
y=yi _ 

The difference equations may be written in the form 

m e n e +c e u e + k e u e ~(\ e P e = i e 
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and 

q e u e =0 

where k e c e , q e , and m e are finite difference “element” matrices, u e is the vector of 
generalized primary variables, ii e is the vector of time derivatives of the generalized 

primary variables, and f e is the finite difference generalized force vector. Assembling the 
element equations yields 

Mu + C(u)u + Ku - QP = F (3.48) 

-Q T u = 0 (3.49) 

where 

nelem T nelem nelem 

K= X 2//a i T9N x=x . ,M= X pN|x=x ; - , C= £ p(Nu e )N|x=x / 

i [ y=y\ J 1 L y=yt -I 1 L y=yi - 

nelem T A 

Q= X 3 IN x=x . , u and u e are vectors that contain all nodal degrees of freedom 

1 . y = yt . 

associated with the primary variables and its time derivative, and 

modes _ , 

F = X I/ 9 *) i x i > , v ; ) + * (- r / ’ yj )J • As in the finite element method, Eq. (3.48) and (3.49) 

1 

can be combined into one system of equations and written in matrix form as 

Oljul rC(u) + K -QlfuljFl 

o oJ[pj + L -Q T o Jjpj joj 

or in a more symbolic form as 

MU + KU = F 


where U={uj U2 U3 P} T . 
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As in the case for the scalar-field formulation, the shape functions for a nine-node 
quadrilateral finite element are used (see Table 2.1). The shape function at point i-lj-l is 
given by 


= |(1 - ^ Xl - 'J) - 1 (l - f 2 )tl - /?) - 1 a - f X> - -7 2 ^ (l - ^ 2 Xl - -7 2 )- 

Similarly, 

N i+\,j-\ = ^( 1 + f Xl- ri)-\{ )(l - *7)"(l + ^ )+ ^t 1 - £ 2 1 1 - n 2 ), 

N M,j+i = ^(1 + iXl + ^)-|(! + 4 1 - n 2 )-j(l -£ 2 )(l + 1 ;) + ^(l -^ 2 )(l - J ] 2 ), 

and 

N i-\,j + i = ^(i - iXi + ^)- 7 (i - £ 2 )(! + n)~ \ (i - Ck ~ ri 2 )+ ^ (i - £ 2 \i - n 2 ), 


Then, for a square element 


dxdy L e V d&T] 

_ a 2 ^ f+u+ i 

dx 2 



( \ 1 * 1 , 

c — n + cn 

>4 2 ’ 2 


> 

>£=0,77=0 



a 2 ^+W— 1 _ 1 ^ 2 7v /+i,y-i 

dxdy ^ e j 2 3£9/7 

a 2 A'-w,i 

dx 2 


1 U 1 . 

C + — 17 + 07 

4 2 2 


d £= 0 , 77=0 


4/7 


(3.50) 


The standard finite difference representation follows by direct substitution of Eqs. (3.50) 
for the cross-derivative terms of the momentum equation along with Eqs. (2.40) for the 
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second-order terms. As noted previously, a single spatial modeling approach (i.e., the 
finite element method or the finite difference method) is used for the single-domain 
formulation. While for multiple domains, homogeneous approaches and heterogeneous 
approaches are available. That is, the same method in each domain (homogeneous 
approach) or different methods in different domains (heterogeneous approach) are 
possible combinations of spatial modeling. 

3.6.3. Multiple-Domain Modeling - Homogeneous Discretization 

These homogeneous approaches make use of a single discretization method 
among all subdomains in which the domain is subdivided. The focus of this work is on 
the finite element and the finite difference methods as the spatial discretization methods. 
For homogeneous domain discretization developed herein, Eq. (3.32) is used to provide 
the mathematical basis for the three-approximation formulation. The generalized element 
equations, for both the finite element and finite difference methods, may be obtained by 
rewriting Eq. (3.32) over an element domain as 


J H)4 1) ) dn i- J Wj a f W - jr/ 2 M 2) dr 1 

f\e t £ t 6 

2 i L h r r 

+ j 4 1 ) ( v <-«? ) )< t i ‘+ |4 2) (v, -»! 2, )dr'' + Ji,-fe (I) +<f , )dr 1 ' 

,e ,e ,e 

r 1 r 1 r 1 

= jA<4 1) -#)4>! 1) dnf+ #4>! 0 dr/ 

“i rf 

+ fA(i>, (2) -vp ) )4 > i 2) dn!+ P, (2 M 2) dlf 


(3.51) 


qs 


1 2 


25 

Note that in the potential energy formulation , the continuity of the secondary variables 
was satisfied through the subsidiary conditions obtained through the minimization of the 
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potential energy. In this weighted residual formulation, the continuity of the secondary 
variables is satisfied in a weighted residual sense and the Lagrange multipliers, Ap ^ and 

X ( , are represented by weight functions in the form of the secondary and primary 
variables, respectively. 

The form of the equations for the finite element and finite difference applications 
differs by the form of the element shape functions and the approximation selected for the 
weight functions, O. The formulation for solid mechanics and fluid mechanics differs by 
the constitutive relations. For the generalized element expansion of subdomain i, the 
independent approximations for the element generalized primary variables, ( i.e ., 
displacements or velocities), interface secondary variables {i.e., tractions or fluxes), the 
weight functions associated with the secondary and primary variables, and the interface 
variables, are, respectively 

u k=^k u e k ’h= R k a k J *-*= R *;^ = T and v = Tu : (3.52) 

Both the solid and fluid mechanics derivations may be developed from Eq. (3.51), given 
the approximations of Eq. (3.52), the appropriate constitutive relation, and the choice of 
weight function. Each derivation is presented in turn in the following work. 

Solid Mechanics- finite element discretization 

Substituting the approximations of Eq. (3.52) into Eq. (3.51) along with the 
constitutive equations and using the Galerkin method in which O = N , yields 
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JB] E]Bj dnf u q + jBiE 2 B 2 daf u e2 + J /?N / N \ dQf ii q + JpNiN 2 dQf 2 ii £ 


- jN^Rjdr 1 " «!- JN 2 R 2 dr r « 2 +| JT^Rj dr 1 " aj + jT T R 2 dr r |a 2 + 


T T T* 


Tn 


jR^Tdr 1 k-l JR/Njdr 1 ^ u +1 jR]Tdr r k-l }R 2 r N 2 dr r u 




Trr 


>Tivt jt^I C 


= JB^Eje® dQf - jB^dflf + jN^pbjdilf + jN/tjdrf' 


+ |BjE 2 £^dQf- jB^d a e 2 + jN 2 pb 2 dQ 2 + jN^drf 


where B^ = 3N^ for /c=l,2 and for the tf h subdomain, the element matrices are 


k e, = / B, T E,B, dQ e k ; m H = j pN T k N k dn< 


k =- jRfcN^dr 1 , 


K k =~ KR k dT l 


kr = fT^dr 1 , 


f ek = JbJe* 4*> dQ|- jBjf«<*>dfl£ + /Nj/3b*dflf+ jN^drf. 
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Solid Mechanics- finite difference discretization 

Substituting the approximations of Eq. (3.52) into Eq. (3.51) along with the 
constitutive equations and using the Dirac delta function as the weight function, 

<t> k =S(x-x i ,y-y i )=S{x i ,y i ), 


3 t e 1 3n 1 

X=Xj 

+ 

3 T E 2 dN 2 

X=Xf 

u e 2 + 

r 

K* 

II 

K 

€ 

i 

u ei + 

I 

K* 

II 

* 

(N 

% 

l 


y=yt _ 



y=yt _ 


L y=yn ] 


L y=vi J 


u 


e 2 


R||x=x z - 

°i- 

jo 

to 

X 

II 

1 

a 2 + 

Jt t R! dr 1 * 

(X| + 

jT T R 2 dr It ’ 

y=yt _ 


L y=yt J 


_r ie 


_r lS 


«2 


+ 


+ 


1 


jR^Tdr 16 

Ui- 

_r lS 

- 

Jr 2 t dr l£ 

Ui- 

_r lS 

- 


dr 1 ' 


u 


jR2N 2 d r 


e l 


e 2 


= 3 T E,4 I) 

+ 3 t E 2E W 

y=y, 

where, for k=\,2 and for the subdomain, the element matrices are 

k e* = J 3 Te /5:3N; 6Q e k - = J pN k \x= Xi dn 


X—Xj 


X=Xj 

X 

II 

H 

+ 

X 

II 

K 

+ 

y=yt 


y=yt 

y=yt 

y=yt 

X=Xj 

-s T 4 2) 

X=Xj 

+ pb 2 \x=xi 

+ h\ x =Xi 

y=yt 


y=yt 

y=yi 

y=yi 


n e 


y=yt 




y=yt 


k P* =- iRlN t dr ie , 


— i x i ’ }’i ) ’ 

k h = jT^dr 1 ', 


(3.55) 


(3.56) 
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and 


l ek 


= 3 T E t e« 


X=Xf 

y=yt 




X=Xj 

y=yt 


+ pb /; \x=Xj + ptk\x=Xi 

y=yi y=yi 


For both the finite element and the finite difference discretization strategies, 
assembling the element equations over the entire domain, enforcing continuity of the 
primary and secondary variables only within each subdomain and assembling the 
contributions along the element edges on the common subdomain boundary, and noting 
that u ej and u e2 , ii e and ii e2 , f ej and f e2 , and aj and <*2 are completely uncoupled. 


yields the system of equations given by 


Mj 

0 

0 

0 

0‘ 


“1 


0 

M 2 

0 

0 

0 


“1 


0 

0 

0 

0 

0 


«I 

► + 

0 

0 

0 

0 

0 


«1 


0 

0 

0 

0 

0 


«2. 



Ki 

0 

0 

K Pi 

0 


0 0 K Sl 

k 2 0 0 

0 0 K h 

0 k? 1 0 

A 1 

K P2 K h 0 


K 

Ki 


s 2 

2 


0 

0 



»1 


fl 


«2 


h 

5 

U I 

> — < 

0 


«1 


0 


P‘2 _ 


0 


( 3 . 57 ) 


or 


r 

O 

O 

§ 


ii 


1 

* 

0 

C/3 

1 


U 


f 

000 

< 

«I 

> + 

O 

O 

< 

«I 

► = < 

0 > 

1 

O 

O 

O 


a 


1 

0 

1 


a 


0 


where K, M, u, and f are the assembled stiffness matrix, mass matrix displacement 


vector, and force vector for the entire structure, and K p , K s , Ki, ui, and a are the 
assembled K p ,, K s ,, Ki., m, anda^ for all interfaces. The assembled stiffness and mass 

matrices, K and M, are block diagonal matrices containing the stiffness and mass 
matrices, K* and M*, of each of the subdomains along its block diagonal. The interface 
“stiffness” matrix thus contains coupling terms that augment the stiffness matrices of the 
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subdomains along the interface. All of the interface “stiffness” terms appear in the 
stiffness matrix with none in the mass matrix. Similar results may be obtained when 
damping is included. As for the scalar-field problem, the three-approximation approach 
for vector-field problems yields systems of equations (see Eqs. (3.57)) of similar form 
and with the same attributes. Again, due to the generalization for the finite difference 
approximations, the system of equations is not necessarily symmetric due to the off- 
diagonal submatrices, K p and K s , nor are they banded or positive definite. Therefore, 
standard Cholesky solvers may not be used, unless full pivoting is performed to obtain 
the solution. The upper diagonal submatrix blocks contain uncoupled stiffness matrices. 
The symmetry of the matrix is determined by the choice of the weight function, O. In 
general, due to the introduction of fictitious nodes for the imposition of boundary 
conditions and loads in the finite difference discretization, the stiffness matrices are not 
symmetric but are positive definite and sparse. The coupling is accomplished through the 
introduction of the coupling terms in the matrices K p and K s for both approaches. 

The number of additional degrees of freedom associated with the interface element is 
generally small in comparison with the total number of degrees of freedom in the 
subdomains. Thus, modeling flexibility is provided at a relatively small computational 
expense. The computational expense in this study may be reduced additionally as the 
efficiency of new solution algorithms for the system of equations in Eq. (3.57) is 
increased. 

While it is convenient to represent the weighted residual form over the domain 
using a single equation, the system of equations, Eq. (3.57) is obtained from the 
individual weighted residual expressions over each of the subdomains and the constraint 
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integrals. The first two matrix equations of the system of equations, Eq. (3.57) are 
derived from the weighted residual statement for subdomain k. That is, 

J oWtrJ^d Q e - \<&ip(b\ k } -v\ k \ d£f - JO^dP^ =0 

n e ’ Q e r s e 

The third matrix equation of the system results from the reciprocity statement of the 
secondary variables. That is, 

+# ) )dr I =o on r 1 . 
r 1 

The fourth and fifth matrix equations result from the continuity requirement for the 
primary variables, which is given by 

° n pI 

r 1 

j4%i -“/^Jdr 1 =0 on r 1 

r 1 

Note that the forms of the coupling element matrices that are not in terms of the 
weight functions are independent of the method of discretization. That is, 

k Pi=- f<N t dr : \ 

r ie 

and 

ki* = l^R^dr 1 " 

T e 

r 1 

are of the same form for the finite element and finite difference discretizations. 


However, since the element shape functions, N*, differ for the two methods, the interface 
matrices, k p ^ , in general, are not identical. 
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Fluid Mechanics- finite element discretization 

Substituting the approximations of Eq. (3.52) into Eq. (3.51) along with the 
constitutive equations and using the Galerkin method in which O = N , yields 


dDj 

+ 


Ug 2 + 

J PN^Npi^dQ* 

Q i 


_°2 


Qj 


J pNj(N 2 u e2 )N 2 dQ^ 

U e2 + 

f 2//B, 1 B lf dQf 

Uq + 

J 2//B 2 B 2f dQ. e 2 

Q 2 


Q l 


Q 2 




jB^N^dQf 

Q i 

P e,- 


JbJ N 2 I dQ^ 


Q 


e 

2 


e 2 




" 


" 



jNjR, dr 1 * 

«l - 

J n Jr 2 dr 1 '' 

« 2 + 

JT 1 Rj dT f 

<X| + 

Jt t r 2 dr 1 * 

_r r * 


J e 

Lr 1 J 


_r ie 


.r 1 * 









jR^Tdr 1 * 

ui - 

JrT Ni dr 1 * 

*ey + 

JrJt dr 1 ' 

«i- 

jR^N 2 dr I<? 

_r ! * 


r l& 

.r 1 * 


S l * 


= jNjpbjdQj + jNj r t 1 drj ?e + |nJ /tb 2 dQ2 + jN 2 t 2 drf 

r f °2 rf * 


where B/. = 9N ^ and B /. =9fN^ for k=\,2 and the elemental matrices are 

k, t = j2//B[B tf dflj; = J /*!>!* ; 

Q| a‘ k 

Ce t = 1 />Nl(N t » et Kdn‘;q et = J B^IcinJ; 


(3.58) 
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k P 4 =- W N * dr 1 * , 
r 1 * 

k *t =- jNjRtdr 1 ', 

r 1 

ki*= jT T R*dr ie , 

i e 


and 


f e k = I N[pMn| + ( N^ t drf. 

Q e r s e 

k r jfc 

In addition to the element equations for momentum, Eq. (3.58), the element 
equations for continuity must be considered. Eq. (3.34) is used to provide the 
mathematical basis for the continuity equation for multiple domains. Using the Galerkin 

A /V 

method, the weight function corresponding to the continuity equation is given by O = N . 
Substituting the approximation for the weight function and the primary variable into Eq. 
(3.34) yields 


jB/NjdO* 


A 



J 


+ 


|BjN 2 dD t 


“2 


M 


^2 


= 0 


Fluid Mechanics- finite difference discretization 

Substituting the approximations of Eq. (3.52) into Eq. (3.51) along with the 
constitutive equations and using the Dirac delta function as the weight function, 

*>k= s k{ x - x i>y-yi)= s k( x i>yi)> 
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/?Nj|x=X,- 

y=yi 


u q + 


pN 2 \x=Xi 

y=yt 


|u e 2 + 


p(Niu q )Ni| x=x . 


y=yt 




+ 


/ 7 ( N 2«e 2 ) N 2|x=x I - 


y=yi 


\u e2 + 


2//a T T3N] 


X=Xf 

y=yt 


u q + 


2jud [ TdN 2 


X=Xf 

y=yt 


e 2 



X=Xi 


a T m 2 

X=Xf 


y=yt _ 



y=yi _ 


e 2 


-[ R 2(W/)]«2 + 

Jt t R! dr 1 * 

«i + 

|T T R 2 dr ie 


_r I<; 


r ie 


«2 


+ 


J Rj" T dr 1 * 

U I“ 

jR^Nidr 1 " 

u q + 

jR^Tdr 1 * 

»I - 

jR 2 N 2 dr ie 

_r J * 


T lC 

.r 1 * 


_r I<? 


u 


= pb l {x i ,y i )+t l (x t ,y i )+ pb 2 {x i ,y i )+t 2 {x i ,y i ) 
where, for k=l,2 and the elemental matrices are 

k C £ — 2/jd T3N^ x =Xf 9 m e^ ~ /^^A:| x=x / > 


y=yt 


*~ e k — k x ~ x i 9 Qet ^ 1^” 


y=yt 


y=yt 


x=x f 9 

y=yt 


k P* =- jR> t dr' 


— ^ k i-^i i) ' 


ki*= jT^^dr 1 


and 


f elr — pb k ( x i > J / / ) ^k^i^yi ) • 


e 2 

(3.59) 
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Considering continuity, and using the Dirac delta function as the weight function, 

<\> k =Sj c (x-x i ,y-y i )=Sj c (x i ,y i ), the continuity equation is given by 

9 Tl Nl| x = X; . « ei + 3 T IN 2 | J c=x / u e 2 = 0 

y=yi J L y=yt . 

For both the finite element and the finite difference discretization strategies, 
assembling the element momentum equations, Eqs. (3.58) and (3.59) over the entire 
domain, enforcing continuity of the primary variable only within each subdomain, and 
noting that u C] and u e2 , u Cl and u e2 , P C| and P e2 , f C| and f e2 , anduj anda 2 are 

completely uncoupled, yields the system of equations given by 

Mi 0 0 0 0 0 Olfui] [Cl 0 0 0 0 0 0][ui 

0 M 2 0 0 0 0 0 u 2 0 C 2 0 0 0 0 0 u 2 

0 0 00000 U! 0 0 00000 uj 

0 0 0 0 0 0 0 < Pi U 0 0 0 0 0 0 0 Pi + 

0 0 0 0 0 0 0 P 2 0 0 00000 P 2 

0 0 00000 «! 0 0 00000 «1 

0 0 0 0 0 0 oJ|a 2 J [o 0 0 0 0 0 oJ|a 2 

0 0 0 Qi 0 0 0 

0 0 0 0 Q 2 0 0 

0 0 0 0 0 0 0 

0 0 0 0 0 0 0 

0 0 0 0 0 0 0 

0 0 0 0 0 0 0 

0 0 0 0 0 0 0 

Ki 0 0 0 0 K Sl 0 

0 K 2 0 0 0 0 K S2 

0 0 0 0 0 0 0 

0 0 0 0 0 0 0 

0 0 0 0 0 K :i Kj 2 

K pi 0 0 0 0 0 

0 K P2 k£ 0 0 0 0 







along with 
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Qi °l{ u i} = {° 

0 Q2_L u 2 J 1° 

or symbolically 

M 0 0 0 ii C(u) + K 0 -Q K s u f 

0 0 0 0 iij 0 0 0 Kj uj 0 

0 0 0 0 P — Q T 0 0 0 | p |o 

0000a Kp 00a 0 

where K, M, C, Q are the assembled coefficient matrices for momentum and continuity, 
u and f are the displacement vector and force vector for the entire structure, and K p , K s , 
Ki, m, and a are the assembled K Pjt , K Sjt , Ku, m, and for all interfaces. 

The first two matrix equations of the system of equations, Eq. (3.60 ) are derived 
from the weighted residual statement for subdomain k. That is, 

J 6£l e - JO,p(^-v^)dD e - JO^dT 5 * =0 

Q e ’ Q e T s e 

The third matrix equation of the system results from the reciprocity statement of the 
secondary variables. That is, 

+^)dr I =0 on r 1 . 

r 1 

The fourth and fifth matrix equations result from the continuity requirement for the 
primary variables, which is given by 

-M^Jdr 1 =0 on r 1 

r 1 

° n pI 

r 1 
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Note that the forms of the coupling element matrices that are not in terms of the 
weight functions are independent of the method of discretization. That is, 

kp* =- (R^dr 1 ', 

r 1 ' 

and 

k lk = jT T R*dr ie 

T e 

r 1 

are of the same form for the finite element and finite difference discretizations. 

However, since the element shape functions, N*, differ for the two methods, the interface 
matrices, k pA _ , in general, are not identical. 

In addition, for the penalty finite element model, the system of equations is of the 

same form as given in Eq. (3.60), except that penalty terms are included rather than the 

pressure terms. The resulting system of equations is given by 

Mj 0 0 0 Oil'll!] rCj 0 0 0 OlTin' 

0 M 2 0 0 0 u 2 0 C 2 0 0 0 u 2 

0 0 0 0 0 \iii U 0 0 0 0 lui > + 

0 0000 «1 0 0 0 0 0 «i 

0 0 o 0 oJ|a 2 J [o 0 0 0 oJ|a 2 

0 0 K s! 0 u i ] [Si 0 0 0 0l[ui 

C 2 0 0 K s 2 u 2 0 S 2 0 0 0 u 2 

0 0 K Ii K I 2 < Ul + 0 0 0 0 0 < 

o k£ 0 o ttl o 0 0 0 0 «1 

; KJ 0 0 a 2 0 0 0 0 0 a 2 

F2 L2 J 




or symbolically 
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M 

0 
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C(u) + K + S 
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K s 
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0 

0 

0 
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«I 

> + 
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Ki 
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U I 
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_ K p 

Kf 

0 


a 
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where K, M, C, S are the assembled coefficient, mass, momentum, and penalty matrices, 
u and f are the displacement vector and force vector for the entire structure, and K p , K s , 
Ki, ui, and a are the assembled K Pjt , K S( ,., K^, m, and a £ for all interfaces. Recall that the 
element penalty matrix for the /c ,h subdomain is given by 

= J ( N lt,( N d*2 dn * 

a? 

3.6.4. Multiple-Domain Modeling - Heterogeneous Discretization 

The multifunctional weighted residual formulation of Eqs. (3.57) and (3.60) are 
used as the mathematical basis for multiple-domain modeling using heterogeneous 
discretization. Considering the two domains upon which this discussion is based, one 
subdomain is discretized using the finite element method, and the other subdomain is 
discretized using the finite difference method. Again, for the finite element development, 
the weight functions for the primary variables, u and P, are taken to be the finite element 

shape functions (i.e., <I>£ = and = N/ f ), and for the finite difference development, 
the weight functions are taken to be the Dirac delta function (i.e., 

(I>£ = Sfc (x — Xj , _y — yi ) = S/ C (x i , T/ )). As expected, the set of element matrices becomes 
a hybrid of the matrices from the finite element method and the finite difference method. 
For completeness, these matrices are repeated here for the finite element and finite 


difference subdomains for solid mechanics as 
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k ei = J bTEjBj dflf and k e2 = 3 T E 2 3N 2 x=x . , 

Qf 7=7/ 

m ei = J pNjNj dfi[ and m e2 = /?N 2 |x=x,- , 

Qf y=y ‘ 

and f e[ = Jb/’Ejs^ dQf - Jb/ o® dQf + jN 1 T y ob 1 dDf+ {N't, drf and 

Qj Qj Qj r 5 e 

fe 2 - 3 T E 2 4 2) -3 T o| ( 2 ^ +/?b 2 |x=Xj +t 2 |x=Xj . 

7=7/ 7=7/ 

For fluid mechanics the element matrices are given by 

k ei = j2//B i r B lf dQj and k e2 = 2//3 T TaN 2 | x=x . , 

Q® 7=7/ 


m e] = J d£2[ and m e2 = /?N 2 |x=x,- , 
flf *=* 

c e, = j P N i 1 (n | )N j dQf andc e2 = />(n 2 u £ , 2 )n 2 | x=x . 

Of 7=7/ 

q ei = jB^NjIdQf and q e2 = 9 T IN 2x=x . , 

Qf 7=7/ 

and f ei = |N 1 r /?b 1 dQ'j ? + jN^tj drf and = pb 2 (x i ,y i ) + t 2 (x i ,y i ). 

Of r/’ 


The coupling matrices at the element level are of the same form for both solid and fluid 
mechanics and these matrices are given by 

k Sl =- jN^Rdr 1 " and k S2 =-R(x / ,y / ) , 


and for the two domains, k= 1,2, 
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k P 4 =-| R t N t dr 1 ', 
r je 

and 

ki*= jT T R^dr ie . 

] e 

r 1 

3.7. COMPUTATIONAL IMPLICATIONS 

The multifunctional modeling approach for the vector field problem has been 
generalized such that it is applicable to solid and fluid mechanics as well as both 
homogeneous and heterogeneous discretization approaches. As such the computational 
implications are presented in this section for the generalized system of equations, Eqs. 
(3.57) and (3.60). Implications specific to a discipline or a discretization approach are 
highlighted, where appropriate. 

The assembled coefficient matrices, K, M, C, and Q, are block diagonal matrices 
containing the matrices, K*, M*, C*, and Q* of each of the subdomains along its block 
diagonal. The interface coupling matrix thus contains terms that augment the coefficient 
matrices of the subdomains along the interface. All of the interface coupling terms 
appear in the coefficient matrix associated with the primary variables with none in the 
matrix associated with the time derivative. Again, due to the generalization for the finite 
difference approximations, the system of equations is not necessarily symmetric due to 
the off-diagonal submatrices, K p and K s , nor are they banded or positive definite. Note 
that, even for a single domain model, the mixed formulation results in a nonpositive 
definite matrix. Therefore, standard Cholesky solvers may not be used, unless full 
pivoting is performed to obtain the solution. The upper diagonal blocks contain 
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uncoupled fluid flow coefficient matrices. The symmetry of the matrix is determined by 
the choice of the weight function, O. For the finite element discretization, the subdomain 
matrices are symmetric. In general, due the imposition of boundary conditions and loads 
in the finite difference discretization, the coefficient matrices, K*, are not symmetric but 
are positive definite and sparse. The coupling is accomplished through the introduction 
of the coupling terms in the matrices K p and K S/ _ for both approaches and each of the 

disciplines discussed herein. 

In addition, due to the generalization for the finite difference approximations, the 
system of equations is not necessarily symmetric due to the off-diagonal submatrices, K p 
and K s . The system unknowns in Eq. (3.57) and (3.60) consist of both primary and 
secondary variables given by the displacements or velocities, u, and the traction 
coefficients, a, respectively. Generally, the coupling matrices, K S/ . , are of the order of 

the length of the interdomain boundary, which results in a marked difference in the 
magnitude of the off-diagonal terms of the system matrix compared to its diagonal terms. 
This characteristic produces an ill-conditioned matrix whose solution can cause 
difficulties for some general-purpose solvers. Hence, the coupling matrix should be 
scaled such that it is of the same order as the subdomain stiffness. 

The load transfer mechanism of the multifunctional approach may be interrogated 
for the vector-field problem by considering the first and second rows of Eqs. (3.57) and 
(3.60) for solid and fluid mechanics, respectively. In either case the matrix equations of 
interest are given for solid mechanics by 

Mjtij +KjUj +K S] ai = fj 
M 2 ii 2 +K 2 u 2 +K S 2 a 2 =f 2 



or for fluid mechanics by 


140 


MjU] + KjU] + Kgjtt[ — fj 

M 2 u 2 +K 2 u 2 +K S2 « 2 =f 2 

These equations can be partitioned such that they correspond only to the primary 
variables, on the interdomain boundary 


K l«l + K Sl «l =0 
K 2 u 2 +K S2 « 2 =0 


(3.62) 


andK ^ denotes stiffness terms related to u /. and there are no forces (including inertial 


forces M / ( ii / ; ) on the interdomain boundary and assuming steady fluid flow (i.e., 

M/.u^ =0). The expressions given by K /. u^. represent the internal fluxes at the 
interdomain boundary, and thus Eq. (3.62) may be written as 

?1 = — K sj «1 and f2=- K s 2 a 2- ( 3 -63) 

For homogeneous discretization using the finite element method, substituting for 
K s ,. from Eq. (3.54) into Eq. (3.63) gives 


f! =- Jn^Rj dr 1 * «j = — jNft! dr 1 " 


(3.64) 


f 2 = Jn 2 r 2 dr ie « 2 =- |Njt 2 dr ie 


(3.64) 


Examining Eqs. (3.64) indicate that the evaluation of the internal forces is consistent with 
the evaluation of equivalent nodal forces in the presence of applied tractions on the 
boundary. In addition, Eq. (3.64) substantiates that the secondary variable along the 
interface is represented by distributed forces for each of the subdomains. 
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For homogeneous discretization using the finite difference method, substituting 
for K^. from Eq. (3.56) into Eq. (3.63) gives 

f] = — R |«] = — 1| (3.65) 

f 2 = - R 2 <*2 = - *2 

Examining Eq. (3.65) indicates that the evaluation of the internal forces is consistent with 
nodal forces evaluated at points in the presence of applied tractions on the boundary. In 
addition, Eq. (3.65) substantiates that the secondary variable along the interface for this 
approach is represented by nodal forces for each of the subdomains. 

For heterogeneous discretization using the combined finite element and finite 
difference methods, substituting for K s from Eq. (3.54) into Eq. (3.56) gives 

fj =- Jn/Rj dr 1 * =— jNftj dr 1 * (3.66) 

,e ,e 

r 1 r 1 

f 2 = -R 2 « 2 = -t 2 

Examining Eq. (3.66)) indicates, for subdomain 1, that the evaluation of the internal 
forces is consistent with the evaluation of equivalent nodal forces in the presence of 
applied tractions on the boundary, while for subdomain 2, the evaluation of the internal 
forces is consistent with nodal forces evaluated at points. This reveals that for this 
multiple domain approach, the secondary variable along the interface for subdomain 1 is 
represented by distributed forces, and for subdomain 2, the secondary variable along the 
interface is represented by nodal or point forces. Thus for this heterogeneous modeling 
approach, it is required to transform the interface secondary variables into equivalent 


quantities. 
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3.8. VERIFICATION TEST CASE 

In this section, the multifunctional methodology for the vector-field problem is 
demonstrated on a verification test case. The application is described and the associated 
results and salient features are discussed. This application is considered a patch test for 
the formulation and verifies the applicability of the method for a configuration for which 
the solutions are known. Finite difference and finite element solutions for single- and 
multiple-domain configurations are presented to provide benchmark solutions for the 
multifunctional approach using heterogeneous discretization. Representative applications 
from the field of engineering science are presented in Chapter V. 

3.8.1. Patch Test 

As in the scalar-field problem, a patch test is used to determine the effectiveness 
of the multifunctional approach applied to a vector-field problem. A cantilevered plate is 
subjected to uniform inplane loading at the free end that yields a constant state of strain. 

In particular, this loading condition provides verification of the finite difference method 
for combinations of displacement and fraction boundary conditions, and the method is 
validated for both the single- and multiple-domain models. 

Problem Statement 

The analysis domain and the boundary conditions are shown in Figure 3.1. The 
normal and tangential fractions are denoted by T„ and T t , respectively, in the figure. This 
configuration has been used in the combined finite difference and finite element analysis 
reported by Dow et al. , and it is used here to provide a point of comparison. The length 
of the plate, L, is 20 in., the width, W, is 8 in., and the thickness, h, is 1 in. The material 
system is described by a Young’s modulus of 30,000 psi and a Poisson’s ratio of 0.3. An 
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applied displacement of 0.3 in. is applied at one end, and the opposite end is fixed. The 
other sides are free. 

For the finite element method, four-node elements are used to discretize the 
domain for all applications. Homogeneous discretizations for single- and multiple- 
domain models for the finite element and finite difference methods are presented. For the 
finite element discretization of a single domain, a finite element mesh of 20 elements and 
4 elements are used in the axial (x-direction) and transverse directions (v-direction) , 
respectively, of the plate. For multiple domains with compatible meshes ( i.e nodal 
coincidence is maintained at the interface), two finite element meshes of 10 elements and 
4 elements are used in the x- andy-directions, respectively. For the finite difference 
discretization of a single domain, a finite difference grid consistent with the finite 
element mesh was used. That is, a grid of 21 grid points and 5 grid points are used in the 
axial (x-direction) and transverse directions (y-di recti on), respectively, of the plate. 
Similarly, for multiple domains with compatible meshes, two finite difference meshes of 
1 1 grid points and 5 grid points are used in the x- andy-directions, respectively. For 
multiple domains with incompatible finite element meshes, one domain is discretized 
with 10 elements in the x-direction and 4 elements in they-direction. While the other 
domain is discretized with 20 elements in the x-direction and 8 elements in the y- 
direction. The multiple-domain discretization is shown in Figure 3.2. The finite 
difference discretization is consistent with the finite element mesh discretization. 
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T„=T,=0 



U=Uo 

T t =0 


Figure 3.1. Analysis Domain and Boundary Conditions of Cantilevered Plate. 


Interface 



Figure 3.2. Multiple- Domain Discretization of Cantilevered Plate. 

Boundary Conditions for Finite Difference Method 

The finite difference method is extensively tested for the single- and multiple- 
domain configurations to assure that the boundary conditions are being applied correctly. 
Generally, for the vector-field problem, a 3x3 or nine-point central difference template is 
used to evaluate the momentum equation, Eq. (3.20). On the boundary of the domain, the 
template introduces fictitious nodes. In reference 43, the fictitious nodes are eliminated 
using traction conditions, T n and T t , and the constitutive equations. When the differential 
equation is evaluated at the comer of the domain boundary (see point i,j in Figure 3.3), a 
fictitious node (point i+lj'+l) is introduced for which there are no additional auxiliary 
equations. Thus, to eliminate the degrees of freedom associated with this fictitious node, 
non-physical higher-order derivatives of the constitutive equations are introduced that 
further complicate the approach. An alternative approach, used herein, is to apply the 
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momentum equation only to the nodes in the interior of the domain, while the differential 
equations representing the traction conditions are applied to the boundary nodes. Special 
forms 44 of the difference equations for grid points at the boundaries are used to avoid the 
use of fictitious nodes. These forms make use of higher-order forward or backward 
difference operators to express the differential forms in order to maintain the same order 
of accuracy as the central difference operator. For multiple-domain spatial modeling, the 
momentum equation is applied to nodes on the subdomain interface boundary. The 
higher-order backward or forward difference operators are used to introduce the unknown 
traction on the interface. This approach yields equations at the interface in terms of the 
unknown tractions at that specified interface node only. If a central difference scheme 
were used for the traction conditions, the equations on the interface would be in terms of 
the unknown tractions at the specified interface node and adjacent interface nodes. In the 
latter case, the resulting equations can not be derived from the generalized 
multifunctional formulation. 


Fictitious Node ■ 


Boundary Node 




Interior Node 


a, 


<? 0 


Tn, T, 


Domain Boundary 




T„,T, 


-6 


Figure 3.3. Central Difference Template Applied at a Corner. 
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Analysis Results 

Several analyses have been performed: (1) two single-domain analyses, one with 
finite element discretization and one with finite difference discretization, respectively, (2) 
two multiple-domain analyses with homogeneous modeling, one with finite element 
discretization in each domain and one with finite difference discretization in each 
domain, and (3) one multiple-domain analysis with heterogeneous modeling with 
combined finite element and finite difference discretizations. All of the analyses yielded 
the exact solution within machine accuracy. Results for the internal forces or stresses 
along the interface for the analysis cases are shown in Table 3.1 . The results are given at 
the locations along the width of the plate normalized by the plate width. 

For the finite element domains, the internal forces, F x and F y , obtained from the 
multiple-domain analyses are normalized by the value of the force obtained from the 
exact solution multiplied by the element length along the edge of the interface. Thus, for 
a consistent load and for the finite elements used in this study, a normalized value of 
unity represents complete agreement with the exact solution at the interior nodes ( i.e ., 
l/8<y/fF<7/8). At the end nodes (i.e.,y/W= 0 andy/lT=l), a normalized value of one half 
represents complete agreement with the exact solution. 

For the finite difference domains, the stresses, o x and T xy , obtained along the 
interface from the multiple-domain analyses are normalized by the value of the normal 
stress obtained from the exact solution. Thus, a normalized value of unity represents 
complete agreement with the exact solution. Values in Table 3.1 for the normalized 
distance along the interface, y/W, annotated with a superscript ‘F’ in parentheses denotes 
results obtained from the most refined subdomain (see Figure 3.2). 
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The single-domain analyses with either finite element discretization or finite 
difference discretization are in excellent agreement with the exact solution. Moreover, 
the interface force and stress results obtained with multiple-domain analyses using 
homogeneous modeling with either finite element discretization or finite difference 
discretization are in excellent agreement with the exact solution. For the heterogeneous 
modeling, the finite difference method was used in the coarsely discretized domain, and 
the finite element method was used in the more refined domain. Note that the stresses are 
used to compare the accuracy of the solution in the finite difference domain, and the 
internal forces are used to compare the accuracy in the finite element domain. The results 
obtained from this heterogeneous modeling approach are in overall good agreement with 
the exact solution. 


Table 3.1. Results of the Multifunctional Approach for the Cantilevered Plate. 


Location 

Analysis Type* 

Along 

SD/FE 

SD/FD 

MD/FE 

MD/FD 

MD/HM 

Interface, y/W 

F x 

Fy 

Ox 

%cv 

F x 

Fy 

Ox 

T xv 

Ox 

F x 

0 . 

0.5 

0.00 

1.00 

0.00 

0.5 

0.00 

1.00 

0.00 

.999 

.499 


1.00 

0.00 

1.00 

0.00 

1.00 

0.00 

1.00 

0.00 

- 

.999 

1/4 

1.00 

0.00 

1.00 

0.00 

1.00 

0.00 

1.00 

0.00 

.999 

.999 

3/8^ 

1.00 

0.00 

1.00 

0.00 

1.00 

0.00 

1.00 

0.00 

- 

.999 

1/2 

1.00 

0.00 

1.00 

0.00 

1.00 

0.00 

1.00 

0.00 

.999 

1.00 

5/W 1 

1.00 

0.00 

1.00 

0.00 

1.00 

0.00 

1.00 

0.00 

- 

1.00 

3/4 

1.00 

0.00 

1.00 

0.00 

1.00 

0.00 

1.00 

0.00 

1.00 

1.00 

7/8^ 

1.00 

0.00 

1.00 

0.00 

1.00 

0.00 

1.00 

0.00 

- 

1.00 

1 

0.5 

0.00 

1.00 

0.00 

0.5 

0.00 

1.00 

0.00 

1.00 

.499 


SD/FE: Single-Domain with Finite Element discretization 

SD/FD: Single-Domain with Finite Difference discretization 

MD/FE: Multiple-Domain with Finite Element discretization 

MD/FD: Multiple-Domain with Finite Difference discretization 

MD/HM: Multiple-Domain with Heterogeneous Modeling (combined finite 

difference and finite element discretizations) 
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CHAPTER IV 

REPRESENTATIVE SCALAR-FIELD APPLICATIONS 

4.1. GENERAL 

In this chapter, the multifunctional methodology is demonstrated on several 
representative scalar-field applications. The governing partial differential equation for 
the scalar-field problem is applicable to a variety of problems in engineering science. A 
sampling of these problems include a torsion problem, a heat conduction problem, and a 
two-dimensional flow problem. The applications are described, and the associated 
multifunctional analysis results and salient features are discussed. Finite difference and 
finite element solutions for single- and multiple-domain configurations are presented to 
provide benchmark solutions for the multifunctional approach using heterogeneous 
spatial discretizations. The finite element models use four-node Lagrange isoparametric 
finite elements, and the finite difference model uses a five-point template to approximate 
the governing differential equation. Stand-alone finite element software is used to 
generate the finite element stiffness matrices. The mathematical computing program 
MATLAB® is used to generate the finite difference matrices and the interface coupling 
matrices and to solve the resulting system of equations. 

4.2, TORSION OF PRISMATIC BAR 

The torsion of a prismatic bar with a rectangular cross-section is used to 
demonstrate the multifunctional capabilities for the Poisson problem. As mentioned in 
Section 2.2. 1 , the torsion problem reduces to the nonhomogeneous partial differential 


equation 
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q. + ?!±=-2Gt> 

tlx 2 dy 2 

in which the stress function, (f> , must be constant along the boundary of the cross-section, 
0 is the angle of twist per unit length of the bar, and G is the shear modulus. The 
configuration of the bar is shown in Figure 4.1, and the analysis domain and the boundary 
conditions, are shown in Figure 4.2. 



Figure 4.1. Prismatic Bar with Rectangular Cross-Section. 

For a solid cross-section, the requirement of a stress-free boundary yields the boundary 
condition, (f> = 0 , on all four bounding surfaces along the bar length. Because of the 
symmetries in the problem, only one quadrant of the rectangular cross-section needs to be 
considered. Moreover, this symmetric model is useful in verifying the application of 
mixed boundary conditions. That is, the application of boundary conditions in terms of 
both primary and secondary variables. The quadrant considered in the symmetric model 
is shown in Figure 4.3. 


The shear stresses in the cross-section are 



150 


T T =-^- 

zx dy ’ zy dx ' 

At the ends of the bar, the first moment integrated over the cross-sectional area must 
equal the twisting moment. This requirement gives 

M t = 2\(p dxdy 

and the twisting moment is related to the angle of twist by 

M t = GJ6 

where J is the torsional constant. 




Figure 4.2. Analysis Domain and Boundary Conditions for Prismatic Bar with 

Rectangular Cross-Section. 


The analytical solution 46 for the stress function is given by 


l = 32G f 2 i -U-ijMV 2 

7T n= 1,3, 5,... « 


cosh(n^y/2a) CQS nm: 
cosh{n7tb ! 2a)\ 2 a 


and by differentiating 
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_ _ 16 GOa £ 1 ( ijn- 1)/2 j 

tt 2 b = i , 3 , 5,...« 2 L 


sinh(n^y/2a) 
cosh(n7tb / 2a) 


cos 


«/a 

2a 


and 


T 


zy 


9^ _ 16(76*7 


dx 


it 


5 4 (_1)( ”" 1)/2 1 

«=1,3,5,... n 


cosh(a#y/2a) 

cosh(«7r6/2a) 


sin 


nm 

2a 


Assuming that b>a, the maximum shearing stress corresponding to the maximum slope, 
is at the middle points (y=0) of the long sides x=±a of the rectangular cross-section. 
Substituting x=a, y=0 and recognizing that 


, 1 1 

1 + — + — + •■• = — 


yields 


‘'max 


16 


n ' 2 n=l,3,5,...« 2 cosh(»^6/2a) 


\G0a . 


In addition, the twisting moment, M t , is given by 


M t = ^-Gd(2af{2b) 


. 192 a ~ 1 , njtb 

1 TT ^ -rtanh— — 

7T b n =\^_ t? 2 a 


V 


J 
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(a) Analysis Domain and Boundary Conditions (b) 6 x 6 Mesh of Grid Points 

Figure 4.3. Analysis Domain, Boundary Conditions and Typical Mesh for One Quadrant 
of Prismatic Bar with Rectangular Cross-Section. 


Spatial Modeling of Prismatic Bar 

Analyses are performed for the case of b=2a (i.e., rectangular cross-section), 
where a and b are dimensions of the cross-section shown in Figure 4.3(a). Three levels 
of grid refinement are used for the spatial modeling, namely meshes of (6 x 6), (1 1 x 11), 
and (21 x 21) grid points, each applied to one quadrant of the domain shown in Figure 
4.3(a). A typical idealization for a (6 x 6) mesh of grid points is shown in Figure 4.3(b). 
Multiple-domain analyses with the spatial modeling of these three levels of grid 
refinement and with coincident nodes along the common subdomain boundary have been 
performed for comparison. For the multiple-domain spatial modeling with non- 
coincident nodes along the common boundary, the mesh discretization of the most 
refined domain is consistent with the discretization used in that same region for the 
single-domain analysis. The mesh in the less refined domain has half the “element” 
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density of that used in the refined domain. This mesh is referred to by the syntax (1 1 x 
1 1 )/(2 1 x 21). The coarse and fine finite element models, shown in Figure 4.4, are used 
in the finite element homogeneous spatial modeling. For the finite difference 
homogeneous modeling and the heterogeneous modeling, a finite difference mesh is used 
that has the same number of grid points as the finite element mesh in the respective 
domain. 



Figure 4.4. Multiple-Domain (llxl 1)/(21 x 21) Idealization. 

Twisting Moment for the Prismatic Bar 

Having found the values of the stress function, <f > , at the grid points in the solution 
domain by the respective spatial discretization approaches, the twisting moment may be 
found by repeated application of the trapezoidal rule for numerical integration. The 
computed twisting moment is then normalized by the analytical solution. The normalized 
twisting moment [M t I ^ t analytical ) Stained using the homogeneous and heterogeneous 

spatial modeling approaches are given in Table 4.1. A value of unity indicates perfect 
agreement with the analytical solution. Results in Table 4.1 indicate that all analyses are 
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in good agreement with the analytical solution. The maximum error in any of the 
computed solutions is less than 6%. The maximum error value for the multiple-domain 
analyses is less than 3% and is observed for the multiple-domain heterogeneous modeling 
analysis (MD/HM) using combined finite difference and finite element discretizations. 
Note that some of the error is intrinsic to the coarse approximation of the integral using 
the trapezoidal rule. The integration error decreases as the mesh refinement is increased. 
A more accurate integration rule such as Simpson’s rule would produce results that are 
more accurate. Independent of the integral approximation, the solution accuracy for each 
of the modeling methods increases as the mesh refinement increases. For the same 
number of nodes or grid points, the finite element discretization yields more accurate 
solutions than the finite difference discretization. The results obtained for the single- 
domain modeling ( e.g ., SD/FE and SD/FD) and the multiple-domain homogeneous 
modeling with coincident nodes along the subdomain boundary are identical or nearly 
identical (see the results for (6 x 6), (1 1 x 11) and (21 x 21) meshes in Table 4. 1). These 
results validate the multifunctional approach for coincident grid points along the 
subdomain boundary. The results obtained for the multiple-domain heterogeneous 
modeling approach with coincident grid points along the subdomain boundary are less 
accurate than corresponding results obtained using homogeneous modeling but are in 
overall good agreement. In addition, with the heterogeneous modeling, the accuracy of 
the twisting moment increases as the mesh refinement increases. With multiple-domain 
modeling using finite element (MD/FE) discretization and with non-coincident nodes, 
the accuracy of the twisting moment is bounded by the accuracy of the less refined (1 1 x 


1 1) and more refined (21 x 21) coincident meshes (see the results for the (llxl 1 )/ (2 1 x 
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21) mesh in Table 4.1). For the multiple-domain finite difference (MD/FD) 
discretization in both domains with non-coincident nodes, the twisting moment is slightly 
less accurate than the results obtained using the (1 1 x 1 1) coincident mesh, which is 
indicative of the error introduced by the finite difference interface constraints along the 
common boundary. For the heterogeneous modeling approach with coincident nodes 
along the interface boundary, the twisting moment is less accurate than the homogeneous 
approach with either finite element modeling or finite difference modeling. These results 
reveal the error introduced in the heterogeneous modeling approach for this problem due 
to the interface constraints. However, recall that the twisting moment is a secondary 
result, and the errors obtained are larger than those obtained for the primary variable, (f > , 
the stress function. For the heterogeneous modeling approach with non-coincident nodes, 
the twisting moment is slightly more accurate than the (11x11) coincident mesh, which 
is indicative of the benefit gained (i.e., more accurate field approximation and interface 
constraint) by the combination of the finite element and finite difference discretizations. 


Table 4. 1 . Normalized Twisting Moment for the Prismatic Bar. 


Analysis 

Type" 

Normalized Twisting Moment, M { /(M t ) analytical 

Mesh Density 

(6 x 6) 

(11 x 11) 

(21 x 21) 

(11 x 1 1 )/(2 1 x 21) 

SD/FE 

0.9871 

0.9944 

0.9976 

- 

SD/FD 

0.9743 

0.9897 

0.9964 

- 

MD/FE 

0.9871 

0.9944 

0.9976 

0.9959 

MD/FD 

0.9746 

0.9898 

0.9964 

0.9834 

MD/HM 

0.9498 

0.9738 

0.9878 

0.9749 


SD/FE: Single-Domain with Finite Element discretization 

SD/FD: Single-Domain with Finite Difference discretization 

MD/FE: Multiple-Domain with Finite Element discretization 

MD/FD: Multiple-Domain with Finite Difference discretization 

MD/HM: Multiple-Domain with Heterogeneous Modeling (combined finite 

difference and finite element discretizations) 
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Maximum Shear Stress for Prismatic Bar 

The maximum shear stress, x max , occurs at x=a and v=0 and is obtained by 
evaluating d<p/dx at that point. For the finite element method, the shear stress may be 
obtained from the element shape functions. However, a more general approximation is 
used herein to compare the finite element and finite difference computations. Generally, 
to determine this partial derivative, d(f>/dx , of the stress function, a smooth curve 
containing the stress function values at the grid points can be assumed to represent the 
function, <j ) . Newton’s interpolation formula 47 , used for fitting such a curve, can be used 
to define the function that is differentiated and evaluated at x=a to give the value of 
maximum shear. However, due to errors introduced in the interpolation for large 
amounts of data, a simple backward-difference approximation with the error of the order 


of Ax was used such that 




d(f> 

V Jx=a,y = 0 


2/^ ^- 2 J 4< ^-lJ + WiJ ) 


where the subscripts, i,j, represent the location of the grid point at which the stress 
function is sampled (i.e., x=a,y=0 in this case) and Ax is the distance between the i"' and 
the i- I th grid point. The values for the maximum shear stress, Xmax, obtained using the 
multifunctional approach with single-domain ( e.g ., SD/FE and SD/FD) and multiple- 
domain analyses are normalized by the analytical solution, and these normalized values 
are given in Table 4.2. A value of unity indicates perfect agreement with the analytical 
solution. The results indicate that all of the analyses are in excellent agreement with the 
analytical solution. The maximum error in any of the computed solutions is less than 2%. 


This maximum error value is obtained for the multiple-domain heterogeneous modeling 
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analysis (MD/HM). In general, the solution accuracy for each of the modeling methods 
increases as the mesh refinement increases. An exception to this characteristic is 
observed for the finite element discretization (see the results for the (11 x 11) and the (21 
x 21) meshes in Table 4.2). In this case, the results are oscillating about the analytical 
solution. For the same number of nodes or grid points, the finite element discretization 
yields more accurate solutions than the finite difference discretization. The results 
obtained for the single-domain modeling and the homogeneous modeling with coincident 
nodes along the subdomain boundary are identical or nearly identical. As in the case for 
the twisting moment, this characteristic indicates that the multifunctional approach does 
not introduce error for the compatible meshes. The results obtained for the multiple- 
domain heterogeneous modeling approach with coincident grid points along the 
subdomain boundary are less accurate than corresponding results obtained using 
homogeneous modeling; however, the results are in overall good agreement. In addition, 
with the heterogeneous modeling, the accuracy of maximum shear stress increases as the 
mesh refinement increases. With multiple-domain modeling using finite element 
discretization and with non-coincident nodes, the accuracy of the twisting moment is 
bounded by the accuracy of the less refined (1 1 xl 1) and more refined (21 x 21) 
coincident meshes (see the results for the (llxl 1 )/(2 1 x 21) mesh in Table 4.2). For the 
finite difference discretization in both domains with non-coincident nodes, the twisting 
moment is slightly less accurate than the (6 x 6) coincident mesh, which is indicative of 
the error introduced by the finite difference interface constraints along the common 
boundary. However, the error for all of the finite difference homogeneous analyses is 
much less than 1%; thus, the difference in the homogeneous modeling is not appreciable. 
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For the heterogeneous modeling approach with non-coincident nodes, the twisting 
moment is slightly less accurate than the (11x11) coincident mesh, which, again, is 
indicative of the benefit gained (i.e., more accurate field approximation and interface 
constraint) by the combination of the finite element and finite difference discretizations. 

Table 4.2. Normalized Maximum Shear for the Prismatic Bar. 


Analysis 

Type* 

Normalized Maximum Shear, analytical 

Mesh Density 

(6 x 6) 

(11 x 11) 

(21 x 2 1 ) 

(llxl 1)/(21 x21) 

SD/FE 

1.009 

0.9997 

0.9993 

- 

SD/FD 

0.9940 

0.9973 

0.9986 

- 

MD/FE 

1.009 

0.9997 

0.9993 

0.9995 

MD/FD 

0.9942 

0.9973 

0.9986 

0.9940 

MD/HM 

0.9842 

0.9904 

0.9948 

0.9902 


SD/FE: Single-Domain with Finite Element discretization 

SD/FD: Single-Domain with Finite Difference discretization 

MD/FE: Multiple- Domain with Finite Element discretization 

MD/FD: Multiple- Domain with Finite Difference discretization 

MD/HM: Multiple- Domain with Heterogeneous Modeling (combined finite 

difference and finite element discretizations) 


4,3. HEAT CONDUCTION PROBLEM 

In this section, the basic equation of heat conduction is described briefly to 
provide a convenient reference for the fundamental concepts and equations governing 
conductive heat transfer. The starting point for heat conduction analysis is Fourier’s law 

48 

given in Cartesian vector form for an isotropic medium 

q = -kVT 

where q is a vector whose components are the heat flow per unit area in the respective 
Cartesian directions, k is the thermal conductivity coefficient that may be a function of 


the temperature, T, and 
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Y7 9 . 9 . 9 

V = — i + — 1 + — k . 

dx dy 9 z 

In an isotropic solid with temperature-dependent thermal conductivity, the law of 
conservation of energy with Fourier’s law yields the thermal energy equation. The law of 
conservation of energy is given by 


9# x , %ly , 9 q 


■ + 


+ - 


+e = pc Hi 


9x 9y 9 z 

V s 

where Q is the internal heat generation rate per unit volume, p is the mass density, c is the 
specific heat, and t is time. For constant thermal properties and steady-state heat transfer, 
the heat conduction problem reduces to a nonhomogeneous partial differential equation 
of the form of Eq. (2. 1 ) and is given by 


-k\ 


( 3 2 T | ^2 

dx 2 + dy 2 


= Q 


In this work, two-dimensional heat conduction in a square plate (see Figure 4.5) is used 
to demonstrate the multifunctional capabilities for thermal analysis. For this problem, the 
time-independent, heat conduction equation is 


"d 2 r | 9 2 7^ 

dx 2 + 9y 2 

subject to the boundary conditions 


= y- inQ = {(x,y): 0<(x,y)<l} 
k 


T = 0 on r p = {lines x = 1 and y = l} 
~)T 

— = 0 on r s = {lines x = 0 and y = 0} 
9 n 
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Figure 4.5. Analysis Domain and Boundary Conditions for the Steady-State Heat 

Conduction in a Square Plate. 

Spatial Modeling of Square Plate 

The spatial discretizations in the analyses were selected to be comparable to those 
reported by Reddy for this problem. Coarse and fine models are used in each of the 
subdomains. The coarse model has a (2 x 3) nodal grid, and the fine model has a (3 x 5) 
nodal grid. The syntax (m x n ) is used to denote spatial modeling with m grid points in 
the x-direction and n grid points in the v-direction. The number of grid points, rather than 
the number of elements, in the coordinate directions are used to describe the mesh 
densities to provide consistency when discussing the finite difference and finite element 
discretizations. Combinations of these mesh densities are used for comparative purposes 
where the letters C and F are used to denote the coarse and fine models, respectively. A 
multiple-domain model with finite element models discretized with a fine (3 x 5) nodal 
grid and a coarse (2 x 3) nodal grid is shown in Figure 4.6. Curves labeled C/C or F/F 
denote multiple-domain coarse or fine models, respectively, with coincident nodes along 
the common subdomain boundary. Multiple-domain analyses with finite element 
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discretization or finite difference discretization are denoted by MD/FE and MD/FD, 
respectively. Similarly, multiple-domain analyses using heterogeneous modeling with 
the combination of finite difference and finite element discretizations are denoted by 
MD/HM. 



Figure 4.6. Homogeneous (3 x 5 ) 1(1 x 3) Idealization. 


Temperature Distribution for Square Plate 

The temperature distribution as a function of the distance along they =6 line is 
shown in Figure 4.7 for the different spatial discretizations and modeling approaches. 
The analytical solution for this problem is given by 



In addition, a 1 -parameter Ritz approximation is given by 
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Results obtained using the multifunctional approach are compared to the 
analytical solution (solid line in the figure) and a Ritz approximation (dashed line in the 
figure). Finite element (see Figure 4.7(a)) and finite difference (see Figure 4.7(b)) 
solutions were obtained using a multiple-domain analysis with homogeneous spatial 
discretization and are in excellent agreement with the analytical solution. These results 
illustrate that the temperature at x=y=0 obtained with a coarse finite difference mesh is 
more accurate than that obtained with a comparable finite element mesh (see curves 
labeled MD/FE-C/C and MD/FD-C/C in Figure 4.7(a) and Figure 4.7(b)). This 
difference decreases as the meshes are refined, although the finite element model 
continues to produce a higher temperature value at x=y=0. The multiple-domain 
analyses with non-coincident nodes produce accurate results even at the subdomain 
common boundaries. The multiple-domain results for heterogeneous spatial 
discretization approaches are shown in Figure 4.7(c) and indicate the effectiveness of the 
multifunctional approach. The fine (3 x 5) nodal grid (see Figure 4.6) is discretized with 
the finite difference method, and the coarse (2 x 3) nodal grid is discretized with the 
finite element method. These results are in overall agreement with the results obtained 
with the homogeneous approaches. The homogeneous and heterogeneous results are 
compared in Figure 4.7(d) for models with non-coincident nodes with a fine model in the 
left domain and a coarse model in the right (see Figure 4.6). These results indicate that 
temperatures obtained with the heterogeneous approach are slightly lower than for the 
homogeneous approach with either finite element or finite difference discretizations. In 
addition, the results, obtained by using the finite difference discretization in one or both 
of the domains, illustrate the slight difference in the temperature at the interface from the 
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left and the right domains. However, note that the uniqueness of the solution along the 
interface boundary is satisfied only in an integral sense and this slight difference does not 
detract from the overall accuracy and effectiveness of the multifunctional approach for 
this Poisson problem. 

An additional analysis has been performed to demonstrate the multifunctional 
capability for an inclined subdomain boundary (boundary not parallel to they- axis). In 
this analysis, multiple-domain modeling with the finite element method is used. The 
finite element model used in the analysis has a (3 x 6) mesh of grid points in the left 
domain and a (2 x 3) mesh of grid points in the right domain as shown in Figure 4.8. The 
results for this multiple-domain finite element analysis are shown in Figure 4.9. These 
results (open squares) are compared to the analytical solution (solid line), the Ritz 
approximation ( dashed line) and the multiple-domain finite element analysis (see Figure 
4.6 for the model discretization) with a subdomain boundary parallel to they-axis (open 
circles). The results indicate the effectiveness of the multifunctional approach for the 


inclined subdomain boundary. 
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(a) Multiple-Domain Finite Element 
(MD/FE) Modeling 


(b) Multiple-Domain Finite Difference 
(MD/FD) Modeling 



(c) Multiple-Domain Heterogeneous (d) Multiple-Domain Homogeneous and 
(MD/HM) Modeling Heterogeneous Modeling 


Figure 4.7. Temperature Distribution Along Insulated Edge of Square Plate. 
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Figure 4.8. Spatial Discretization for Inclined Interface for Square Plate. 



Figure 4.9. Temperature Distribution Along Insulated Edge of Square Plate with Inclined 

Interface. 


4.4. POTENTIAL FLOW PROBLEM 


A two-dimensional fluid flow problem is used to demonstrate the multifunctional 
capabilities for a fluid mechanics problem. As shown in Section 2.2.2, the equation 
governing irrotational fluid flow reduces to the Laplace equation 
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d 2 u d 2 u 

— T + —r = 0 

dx 2 dy 2 

where u can be either the stream function, \| /, or the velocity potential, (ft . In this work, 
the two-dimensional, steady, inviscid flow between two infinite plates is considered. A 
rigid, infinite cylinder or radius, R, with an axis at a right angle to the flow is assumed to 
be in the passageway between the plates as shown in Figure 4. 1 0. Far upstream from the 
cylinder there is a uniform flow field with a velocity of Vo. Because of the symmetries in 
this problem, only one quadrant of the domain is considered. The analysis domain and 
the boundary conditions on the velocity potential, <f > , are shown in Figure 4. 1 1 . 



Figure 4.10. Domain of Flow Around Cylinder. 
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Figure 4.11. Analysis Domain of Flow Around Cylinder. 


The finite element models used in this problem are shown in Figure 4.12. A reference 
solution is obtained using the finite element model shown in Figure 4.12(a). The local 
and global finite element models used in the homogeneous and heterogeneous spatial 
modeling approaches are shown in Figure 4. 12(b). For the heterogeneous modeling, a 
finite difference mesh is used in the coarsely refined domain that has the same number of 
grid points as the finite element mesh used in the same domain. This discretization 
strategy illustrates the use of the finite element method to represent the complex 
geometry around the cylinder and the use of the finite difference method away from the 
curved boundary where it is most suitable. 
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(a) Reference Model 


(b) Multiple Domain Model 


Figure 4.12. Spatial Discretization for One Quadrant of Domain of Flow Around 

Cylinder. 

Contour plots for the velocity potential, the horizontal velocity component and the 
transverse velocity component are shown in Figure 4. 13, Figure 4. 14, and Figure 4.15, 
respectively. In each of these figures, the results using the multifunctional approach are 
compared to results obtained from the single-domain analysis using the reference model 
(see Figure 4. 12(a)). As shown in the figures, the velocity potential and the velocity 
components obtained using the multifunctional approach are in excellent agreement with 
the reference solution. In the multiple-domain analyses, the slight discontinuity in the 
horizontal and transverse velocity components at the interface (see Figure 4. 14(b) and 
Figure 4. 15(b)) is due in part to the difference in the computation of the velocity across 
the interface. Unlike in the single-domain analysis (i.e., reference solution), in the 
multiple-domain analyses, the velocities are not averaged across the interface. 
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(a) Single-Domain Model 


(b) Multiple-Domain Model 


Figure 4.13. Contour Plot of Velocity Potential for Flow Around Cylinder. 




(a) Single-Domain Model 


(b) Multiple-Domain Model 


Figure 4.14. Contour Plot of Horizontal Velocity Component for Flow Around Cylinder. 
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(a) Single-Domain Model (b) Multiple-Domain Model 

Figure 4. 15. Contour Plot of Transverse Velocity Component for Flow Around Cylinder. 


The analytical potential solution for the tangential velocity around a cylinder in an 
infinite domain, valid on the cylinder surface, is given by 

f ->2 \ 


«/ = v 0 


1 R ‘ 
1+ — 

r 


sin# 


V J 

where the angle, 0, radial distance, r, and the tangential velocity, u t . can be computed 
from the relations 


6 = tan 


-1 


y 


.a-Xj 


r = 


( a - x) 2 + y‘ 


1/2 


U( = u x sin 6 + u y cos 0. 


The tangential velocity as a function of the angular distance along the cylinder surface is 
shown in Figure 4. 1 6. Results are shown for the tangential velocity around a cylinder in 
an infinite domain for which an analytical solution is known and in a finite domain for 


which a reference solution is obtained using a refined single-domain finite element 
model. For the infinite domain configuration, the plate length to cylinder radius ratio 
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2a/ R , and the plate width to cylinder radius ratio, 2b/ R, are 40 and 20, respectively, and 
the domain can be considered as infinite. That is, the cylinder radius, R , is very small 
compared to the length, 2a, and the width, 2b. For the finite domain configuration, the 
plate length to cylinder radius ratio, 2a/ R , and the plate width to cylinder radius ratio, 

2b/ R , are 4 and 2, respectively, and the domain is considered to be finite. The tangential 
velocity obtained for the multifunctional approach is in overall good agreement with the 
analytical solution for the infinite domain and with the reference solution ( i.e ., single- 
domain analysis) for the finite domain. Results obtained with homogeneous multiple- 
domain analyses with finite element discretization in each domain are denoted by open 
circles in the figure. Results obtained with heterogeneous multiple-domain analyses with 
combined finite difference and finite element discretization are denoted by open squares 
in the figure. The tangential velocity obtained with the homogeneous modeling approach 
is in excellent agreement with the analytical and reference solutions for the infinite and 
finite domain configurations. The tangential velocity obtained with the heterogeneous 
modeling approach is more accurate for the infinite domain configuration than for the 
finite domain configuration. This characteristic is indicative of the performance of the 
finite difference approach, for this problem, in a gradient region. 
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Angle, 0 


Figure 4. 16. Tangential Velocity for Flow Around Cylinder. 

4,5. SUMMARY 

In this chapter, the multifunctional methodology has been described and 
demonstrated for a variety of problems in engineering science. These selected problems 
included second-order problems of solid mechanics, heat transfer, and fluid mechanics 
that can be formulated in terms of one dependent variable. The governing equation in 
each case is either the Laplace or the Poisson equation. The analyses performed have 
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demonstrated the effectiveness and accuracy of the solutions obtained for the respective 
problems. In all cases, the results obtained using the multifunctional methodology were 
in overall good agreement with the reported analytical or reference solution. In the next 
chapter, the multifunctional methodology is demonstrated for problems whose motion is 
described by coupled partial differential equations expressed in terms of two dependent 
variables — vector-field problems. 
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CHAPTER V 

REPRESENTATIVE VECTOR-FIELD APPLICATIONS AND 

EXTENSIONS 

5.1. GENERAL 

In this chapter, the multifunctional methodology is demonstrated on two 
representative vector-field applications. The applications are described and the 
associated results and salient features are discussed. The applications include a plane 
stress problem and a plane flow problem. Finite difference and finite element solutions 
for single- and multiple-domain configurations are presented to validate the 
multifunctional approach using heterogeneous discretization. The finite element models 
use four-node Lagrange isoparametric finite elements, and the finite difference model 
uses a nine-point template to approximate the governing differential equation. Stand- 
alone finite element software is used to generate the finite element stiffness matrices. 
The mathematical computing program MATLAB® is used to generate the finite 
difference matrices and the interface coupling matrices and to solve the resulting system 
of equations. In addition, extensions to multiple discipline analyses are discussed. 

5.2, PLANE STRESS PROBLEM 

A rectangular plate of uniform thickness subjected to a uniform tensile load and 
with a central circular cutout (shown in Figure 5.1) is an ideal example problem with 
which to verify the multifunctional approach. The example problem has a variety of 
practical applications (i.e., rivet holes, aircraft door and window openings, etc.), and an 
exact solution is available 46 . The plate has been used by many researchers to verify 
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aspects of proposed computational methodologies. For example, the plate problem has 
been used by Ransom to verify global/local analysis technology, by Aminpour et al. to 
verify multiple-domain homogeneous modeling using the finite element method, and by 
Rose 10 to verify an adaptive geometry generator used with a multiple-domain finite 
element model. The plate configuration is such that the state of stress is represented by 
the condition of plane stress or plane strain. The membrane displacements, u and v, in 
the axial (x-direction) and transverse (y-direction) directions, respectively, represent the 
plate configuration in plane stress and plane strain. 

Two configurations of this problem have been studied: an infinite plate and a 
finite-width plate. The infinite plate configuration has a central cutout that is very small 
relative to the length and width of the plate, and the exact solution for this problem was 
obtained by Timoshenko 46 . The stress distribution in the neighborhood of the cutout 
exhibits a stress concentration, but from Saint-Venant’s principle, the stress distribution 
is essentially uniform at distances that are large compared with the radius of the cutout. 
The finite-width plate configuration has a larger central cutout relative to the length and 
width, and the stress distribution away from the cutout is not uniform. The finite-width 
plate with a central circular cutout has been discussed by Howland 49 and Peterson 50 . 

For the infinite plate configuration, herein, the length to radius ratio, 2a/ R , and 
the width to radius ratio, 2b/ R , are 40 and 20, respectively, and the plate can be 
considered as infinite. That is, the cutout radius, R, is very small compared to the length, 

H 

2a, and the width, 2b. The material system is aluminum with a Young’s modulus of 10 
psi, and a Poisson’s ratio of 0.3, and the thickness of the plate, h, is 0. 1 in. A uniform 
running load, (N x )o, is applied to each of its ends, and the other sides are free. The plate 
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example problem is used to verify the multifunctional approach for both homogeneous 
and heterogeneous spatial modeling. Because of the symmetry that exists, one quadrant 
of the domain (see Figure 5.2) is modeled. In addition, boundary conditions are shown in 
Figure 5.2 where T n and T t denote normal and tangential tractions, respectively. For the 
multiple-domain analysis, a refined model is used in the near-field subdomain (i.e., the 
local region near the cutout), and a coarse, less-refined model is used in the remainder of 
the domain. A single-domain analysis using a finite element model that has the same 
number of nodes and elements in the near-field region as the multiple-domain model is 
used to obtain a reference solution with which to compare the solution obtained with the 
multifunctional approach. The single-domain model and the multiple-domain model 
(used in the homogeneous spatial modeling) are shown in Figure 5.3. For the 
homogeneous modeling, a finite element (FE) mesh is used in each region. For the 
heterogeneous modeling, a finite difference (FD) mesh is used in the far-field region that 
has the same number of grid points as the finite element mesh in that region. A finite 
element mesh is used in the region near the cutout. 



W 0 =( a x)o fl 


Figure 5.1. Domain of Plate with Central Circular Cutout. 
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T„=T t = 0 



Figure 5.2. Geometric Configuration for One Quadrant of Plate with Central Circular 

Cutout. 



(a) Single-Domain Model 


(b) Multiple-Domain Model 


Figure 5.3. Finite Element Models for One Quadrant of Infinite Plate with Central 

Circular Cutout. 


The exact elasticity solution 46 for an infinite plate with a circular cutout loaded in 
tension indicates that the stress concentration factor, K u is 3.0 at the edge of the cutout 
and is given by 

,= W ,)0 ■ 

The stress concentration factor is defined as the ratio of the maximum stress resultant, 
(N x )max, to the uniform far-field stress resultant, (N x )o- Stress concentration factors 
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obtained using the multifunctional approach with homogeneous and heterogeneous 
modeling are 3.08 and 3.10, respectively, which is within 2.7% and 3.3% of the elasticity 
solution. The stress distributions of the hoop stress resultant (N x )e along the midwidth, 
0=7i, (denoted as line AB in Figure 5.2) and midlength, 9=7r/2, (denoted as line CD in 
Figure 5.2) normalized by the far-field stress resultant (N x )o, are shown in 
Figure 5.4 as a function of the distance from the plate center normalized by cutout radius, 
R. The elasticity solution for the stress distribution is given by 

v Jf, R 2 ) i 37? 4 1 ^ 

)e -t(n x )o 1+_ 2 — +— cos2 ^ 

/ r r 

LV J \ J 

and is shown by the solid lines in the figure. The stress distributions obtained from the 

multifunctional analyses using homogeneous modeling are indicated by the open circles 

in the figure. The stress distributions obtained from the multifunctional analyses using 

heterogeneous modeling are indicated by the open squares in the figure. Excellent 

correlation is observed for all analyses. 

Contour plots of the magnitude, 8 , of the displacement vector ( i.e 

I 2 2 

8 = yu +v ) superimposed on the deformed shape and the longitudinal stress 
resultant, N x , are shown in Figure 5.5 and Figure 5.6, respectively. The multiple-domain 
analysis results are shown for homogeneous modeling using finite element discretization 
in each of the subdomains. To aid visual comparison, the deformation has been 
magnified by 10% of the maximum domain dimension. The displacement contour plots 
reveal the nearly linear variation along the plate length in the far-field region of the plate 
with only local changes near the cutout. The stress resultant contour plots reveal the 
uniform stress state away from the cutout and the peak stress in the neighborhood of the 
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cutout. While not shown, the results for the multiple-domain heterogeneous modeling 
approach are nearly identical to those shown in Figure 5.5 and Figure 5.6, and thus have 
not been included. These contour plots illustrate further the excellent correlation among 
the multifunctional approach using homogeneous and heterogeneous modeling and the 
single-domain solution. 



Figure 5.4. Longitudinal Stress Distribution along Midwidth and Midlength for Infinite 

Plate with Central Circular Cutout. 
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(a) Single-Domain Model (b) Multiple-Domain Model 

Figure 5.5. Displacement Magnitude Distribution for Infinite Plate with Central Cutout. 



(a) Single-Domain Model (b) Multiple-Domain Model 

Figure 5.6. Longitudinal Stress Resultant Distribution for Infinite Plate with Central 

Cutout. 


While the infinite plate analyzed, herein, is an excellent test of the multifunctional 
approach, gradients in the deformation and the stress resultants, as indicated in Figure 5.5 
and Figure 5.6, are well away from the subdomain interface boundary. Thus, to assess 
the accuracy of the approach when the subdomain interface is within a high gradient 
region, a second configuration is analyzed. 

In the finite-width plate configuration, the length to radius ratio, 2a/ R , and the 
width to radius ratio, 2b/ R , are 4 and 2, respectively, and the plate is considered to be 
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finite. The aluminum material system and the thickness that was used for the infinite 
plate is used here for the finite-width plate. The finite-width effects on the stress 
concentration factor for isotropic plates with cutouts have been reported by Peterson 50 . 

By including finite-width effects, the stress concentration factor is reduced from the value 
of three for an infinite plate. The stress concentration factor should be applied to the 
nominal stresses, which are based on the net cross-sectional area associated with the load 
application. For the case of a finite-width plate with a cutout, the net cross-sectional area 
corresponds to 

A net ={ 2b-2R 0 )h = 2bh(\-lf 

{ b 

where h is the plate thickness, and the nominal longitudinal stress for an uniaxial load, P, 
can be expressed as 

p 

(cr v ) = and (iV v ) = (cr v ) h . 

' x / nom a X /nom ' x /nom 

A net 

The geometry definition for the finite-width plate, herein, gives a stress concentration 
factor of 2. 1 6 reported by Peterson. 

Multiple-domain homogeneous and heterogeneous modeling approaches are used 
for the finite-width plate. A refined model is used in the near-field domain, and a less- 
refined model is used in the far-field domain. The single-domain model and the multiple- 
domain model are shown in Figure 5.7. In the multiple-domain homogeneous modeling 
approach, finite element (FE) discretization is used in each domain. In the multiple- 
domain heterogeneous modeling approach, finite difference (FD) discretization is used in 
the far-field domain, and finite element (FE) discretization is used in the near-field 


domain around the cutout. 
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(a) Single-Domain Model (b) Multiple-Domain Model 

Figure 5.7. Finite Element Models for One Quadrant of Finite-Width Plate with Central 

Cutout. 


Stress concentration factors obtained using the multifunctional approach with 
homogeneous (multiple-domain FE/FE) and heterogeneous (multiple-domain FD/FE) 
modeling are 2.19 and 2.73, respectively. These factors are higher by 1.4% and 26.4%, 
respectively, than the values given in Peterson 46 . Note that the solution obtained using 
the heterogeneous modeling approach with finite difference and finite element 
discretizations is nearly 30% in error. This error is likely due to the inaccuracy of the 
finite difference method in the high gradient region and to the constraint conditions along 
the interface. 

To delineate this error, additional heterogeneous analyses are performed using 
finite difference domains with grid spacing in the transverse direction of one half (i. e . , 
9x9 mesh of grid points) and one fourth (i.e., 17x17 mesh of grid points) the grid 
spacing in the initial finite difference domain (i.e., 5x5 mesh of grid points) (see Figure 
5.7(b)). The stress concentration factors obtained with these more refined finite 
difference discretizations are 2.42 and 2.31, which are within 12% and 7% of the 
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Peterson’s solution 46 . The stress distributions of the hoop stress resultant (N x )o 
midlength, 0=7t/2, (denoted as line CD in Figure 5.2) normalized by the nominal stress 
resultant (AQnom, are shown in Figure 5.8 as a function of the distance from the plate 
center normalized by cutout radius, R. The analytical solution reported by Howland 49 is 
denoted by the thick solid line. The stress distribution obtained using 5x5, 9x9 and 
17x17 mesh of grid points are denoted by the short dashed line, the thin solid line, and 
the dashed and dotted line, respectively. The results shown in Figure 5.8 indicate that the 
error decreases as the finite difference grid density increases, and the error decreases 
away from the edge of the cutout. 

The stress distributions of the hoop stress resultant (N x )q along the midwidth, 0=7r, 
(denoted as line AB in Figure 5.2) and midlength, 0=7i/2, (denoted as line CD in Figure 
5.2) normalized by the nominal stress resultant (N x )m,m, are shown in Figure 5.8 as a 
function of the distance from the plate center normalized by cutout radius, R. The 
analytical solution reported by Howland 49 is denoted by the solid lines. This analytical 
solution is valid for distances, r, away from the cutout of less than the plate half-width, b. 
Thus, for this configuration the solution along the midwidth is valid only for r < 2 R. The 
stress distribution for the multifunctional analysis using homogeneous modeling with 
finite element discretization in each of the domains is denoted by the open circles in the 
figure. The stress distribution for the multifunctional analysis using heterogeneous 
modeling with combined finite difference and finite element discretizations is denoted by 
the open squares in the figure. For the heterogeneous modeling approach, the 
distribution is given for the most refined finite difference discretization (i.e., 17x17 mesh 
of grid points). The stress distributions obtained with the multifunctional approach using 
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homogeneous and heterogeneous discretization are in excellent agreement with the 
reported solution. 

Contour plots of the magnitude of the displacement vector superimposed on the 
deformed shape and the longitudinal stress resultant, N*, are shown in Figure 5.10 and 
Figure 5.11, respectively. Results for the multiple-domain homogeneous modeling 
approach using finite element discretization in each of the subdomains are shown in the 
figures. While not shown, the results for the multiple-domain heterogeneous modeling 
approach are nearly identical to those shown in Figure 5.10 and Figure 5.11, and thus 
have not been included. Note that the deformation has been magnified by 10% of the 
maximum domain dimension. The displacement contour plots reveal a deviation from 
the nearly linear variation observed in the far-field region of the infinite plate, and the 
deformation at the cutout is more pronounced. The contour plots illustrate further the 
excellent correlation of the deformation (primary variable) patterns predicted using the 
multifunctional approach with the single-domain solution even with the interface 
boundary domain in a high-gradient region. The stress resultant (secondary variable) 
patterns predicted using the multifunctional approach are also in excellent agreement. 
The slight discontinuity in the stress resultant at the subdomain boundary ( i.e interface) 
is due to the derivation of the nodal stress resultant values from the element quantities. 
The stress resultants are recovered at the finite element nodes by extrapolating the 
stresses at the integration points to the nodes. A single nodal value of the stress resultant 
is obtained by averaging the stress resultants of the adjacent elements. In the multiple- 
domain analyses, the stress is not averaged across the subdomain boundary; thus, any 
gradient across the interface is not considered. 
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Figure 5.8. Convergence of Longitudinal Stress Distribution along Midlength for Finite- 
Width Plate with Central Circular Cutout. 



Figure 5.9. Longitudinal Stress Distribution along Midwidth and Midlength for 
Finite- Width Plate with Central Circular Cutout. 
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(a) Single-Domain Model (b) Multiple-Domain Model 

Figure 5.10. Displacement Magnitude Distribution for Finite-Width Plate with Central 

Circular Cutout. 



(a) Single-Domain Model (b) Multiple-Domain Model 

Figure 5.1 1. Longitudinal Stress Resultant Distribution for Finite-Width Plate with 

Central Circular Cutout. 



5.3. PLANE FLOW PROBLEM 

The flow of a viscous incompressible material squeezed between two long 
parallel plates 41 is considered to illustrate the applicability and performance of the 
multifunctional approach to a representative vector-field problem in fluid mechanics. 
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The geometric configuration and the associated boundary conditions of the problem are 
indicated in Figure 5.12. 

A state of plane flow exists when the length of the bounding plates is very large 
compared to the width of and distance between the plates. Assuming the conditions of 
plane flow, the velocity and pressure fields are determined for a fixed distance between 
the plates. The plates are moving toward each other with a velocity, Vo, and the width of 
and distance between the two plates is given by 2 a and 2b, respectively. For this 
configuration, the ratio of the plate width and the distance between the plates, 2a/2b, is 3. 
Due to the double symmetry present in the problem, one quadrant of the domain was 
analyzed. The viscosity, p, of the fluid is 1 lb-sec/in . The penalty finite element 
model is used in the analysis. The penalty function formulation (see Eq. (3.61)) 
involves treating the continuity equation as a constraint among velocity components. A 
10 x 6 nonuniform mesh (10 elements in the x-direction and 6 elements in they- 
direction) of four-node bilinear elements is used for the single-domain analysis (i.e., 
reference model in Figure 5.13(a)). The nonuniform mesh, with smaller elements near 
the free surface at x=a, is used to delineate the singularity in the shear stress at the point, 
x=a, y=b. This singularity and the associated necessity for nonuniform mesh refinement 
make this problem ideal for demonstrating the multifunctional approach with detailed 
local modeling. The finite element models for the single- and multiple-domain analyses 
are shown in Figure 5.13(a) and Figure 5.13(b), respectively. In the multiple-domain 
analysis, homogeneous spatial modeling with finite element discretization is used. In this 
analysis, more elements are used in the region near x=a, y=b than in the single-domain 
analysis (see Figure 5.13). This local modeling yields a more complex configuration of 
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the subdomain common boundary. That is, the interface between the subdomains 
consists of two non-collateral segments. 


y 



T n = 0 
T,= 0 


-► 

x 


bounding 

plate 


Figure 5.12. Geometric Configuration for Fluid Squeezed Between Parallel Plates. 



(a) Single-Domain Model 
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x=2a/3 x=5a/6 x=a 
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(b) Multiple-Domain Model 


Figure 5.13. Finite Element Models for Fluid Squeezed Between Two Parallel Plates. 


An approximate analytical solution to this two-dimensional problem is provided 


by Nadai 51 and is given by 
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Note that this approximate solution does not satisfy the traction-free conditions ( T n -a x =0 
and T t - t xy =0) on the free edge (i. e . , x=a). Likewise, these traction- free conditions are 
not imposed in the finite element analysis; thus, the conditions are not identically 
satisfied. The horizontal velocity, u, as a function ofy, at three representative locations, 
x=2a/3, x=5a/6 (along the vertical interface), andx=a, is shown in Figure 5.14(a), Figure 
5. 14(b), and Figure 5. 14(c), respectively. The analytical solution of Nadai 51 is 
represented by the solid line in the figure. Finite element solutions obtained using a 
single-domain spatial discretization are represented by the dashed lines in the figure. The 
multiple-domain results for the homogeneous spatial modeling approach using finite 
element discretization in each of the subdomains are also shown in the figures, and these 
results are represented by the open circles. The results for the horizontal velocity 
component obtained from the single- and multiple-domain analyses are in excellent 


V = M\ 


du dv 
dy dx 



190 


agreement with each other, and the results are in overall good agreement with the 
analytical solution. 

The pressure, P, as a function of x, near the centerline for the flow (i.e., y = b/1 6 - 
the centroids of the first row of finite elements in Figure 5.13), is shown in Figure 5.15. 
The analytical solution is denoted by the solid line. The solutions obtained from the 
single- and multiple-domain analyses are denoted by the dashed line and open circles, 
respectively. The results obtained from the multiple-domain analysis are in excellent 
agreement with those from the single-domain analysis. These finite element results are 
also consistent with the results published in the literature . However, the finite element 
models predict a higher pressure in the center of the flow field (i.e., x=0) than predicted 
by the analytical solution. 

While the velocity components and pressure field characterize the flow through 
the plates, the shear stress distribution illustrates the significance of using a graded 
single-domain mesh and a locally-refined multiple-domain mesh. The shear stress, %, as 
a function of x, near the upper bounding plate (i.e., y = 15b/l 6 - the centroids of the last 
row of finite elements in Figure 5.13), is computed at the center of the finite elements and 
is shown in Figure 5.16. Again, the single-domain (dashed line in the figure) and 
multiple-domain (open circles in the figure) results are in excellent agreement with the 
approximate solution of Nadai 51 (solid line in the figure) away from the free-edge. In 
addition, because of the local refinement at the free edge, the multiple-domain results for 
x>5a/6 correspond to the shear stress located at y = 3 lb/32 (the centroids of the last row 
of elements in the refined region). These results illustrate the better representation of the 
gradient in the shear stress at the free edge than either the single-domain analysis or the 
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analytical solution. The approximate nature of the analytical solution is highlighted by 
these results since the solution given does not delineate the gradient on the boundary. 

For completeness, the longitudinal, o Y , and transverse stress, a y , distributions are 
shown at y = 15b/16 and v = 31b/32, respectively, in Figure 5.17 and Figure 5.18. In 
general, the stresses predicted by the single- and multiple-domain finite element analyses 
have a larger value than those obtained by the analytical solution. However, the 
analytical solution is an approximate solution, and the finite element solutions predict the 
same overall trends in the stress distributions. The longitudinal stress distribution, o x , 
reveals the oscillatory nature of the finite element solution at the free-edge. The 
wavelength of the oscillations decreases as the mesh is increased in the local region at the 
free edge as indicated by the results from the multiple-domain analysis. In addition, the 
value of the peak stress at the free edge increases as the finite element mesh is refined. 

Overall, the results obtained with the multifunctional discretization approach are 
in excellent agreement with the single-domain analysis results and with the analytical 
solution given in the literature. These successful comparisons indicate the effectiveness 
of the method and its applicability to the vector-field problem, specifically that of the 


fluid flow problem. 
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(a) Velocity, u, at x=2a/3 



(b) Velocity, u, atx=5a/6 



(c) Velocity, u, at x=a 


Figure 5.14. Horizontal Velocity for the Flow Between Two Parallel Plates. 
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Figure 5.15. Pressure Distribution Near Centerline for the Flow Between Two Parallel 

Plates. 



x/a 


Figure 5.16. Shear Stress Distribution Near Plate Boundary for the Flow Between Two 

Parallel Plates. 
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x/a 


Figure 5.17. Longitudinal Stress Distribution Near Plate Boundary for the Flow Between 

Two Parallel Plates. 


a r (x/a, 3 lb/32) 



Figure 5.18. Transverse Stress Distribution Near Centerline for the Flow Between Two 

Parallel Plates. 
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5.4. EXTENSIONS TO MULTIPLE DISCIPLINES 

In the present work, the multifunctional capability has been demonstrated on 
scalar- and vector-field problems applicable to the general field of engineering science 
and mechanics. While the demonstrations have illustrated the capability within different 
disciplines (i.e., solid mechanics, fluid mechanics, and heat transfer), the method’s use 
has not been demonstrated for multidisciplinary analysis. Extensions to simultaneous 
multiple disciplines are discussed here. 

The term multidisciplinary or coupled systems refers to two or more systems that 
interact with each other, with the independent solution of any one system being 
impossible without simultaneous solution of the others . In general, coupled systems 
and formulations, such as the multifunctional methodology presented in this work, are 
those applicable to multiple domains and dependent variables which usually describe 
different physical phenomena, and in which (1) neither domain can be analyzed 
independently; and (2) neither set of dependent variables can be explicitly eliminated at 
the differential equation level. The class of coupling problems that are the focus of this 
work can be categorized by coupling that occurs on domain interfaces via the boundary 
conditions imposed on that interface. Generally, the domains describe different physical 
situations, but it is possible to consider coupling between domains that are physically 
similar in which different discretization strategies have been used. Fluid-structure and 
thermal-structure interaction problems are typical examples that involve different 
disciplines in different but adjacent domains. Structure-structure or fluid- fluid interaction 
problems are examples where the interface divides arbitrarily chosen regions in which 
different mathematical approximations and/or spatial discretization procedures are used. 
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Single discipline interaction problems have been demonstrated extensively in this work. 
The extension of the multifunctional approach to multiple disciplines is illustrated using 
the fluid-structure interaction problem. 

Different methodologies have been developed for the computational analysis of 
the fluid-structure interaction problem, and different terminology has been used to 
describe the extent to which the disciplines are coupled. In this work, two classes of 
coupling are outlined; namely, fully coupled and loosely coupled methods. Fully coupled 
methods reformulate the governing equations so both the fluid and structural equations 
are combined into one set of equations, coupling the solutions only at the boundary 
interfaces between the fluid and the structure"’ 6 . These new governing equations are 
solved and integrated in time simultaneously. Loosely coupled methods make use of 
independent computational fluid dynamic (CFD) and computational structural mechanics 
(CSM) software modules. The coupling is accounted for by the exchange of data at the 
interface between the fluid and the structure. This coupling approach takes full 
advantage of the numerical procedures of individual disciplines such as finite difference 
approximations for fluids and finite element approximations for structures. In addition, 
software development efforts are simplified and software modularity is preserved. An 
alternate to the coupling approaches is to solve both the structures and fluids problems in 
a single computational domain. The major disadvantage of this methodology is the ill- 
conditioned matrices associated with the two physical domains. A secondary 
disadvantage is the inability to use existing CFD codes because they do not account for 
the interaction with the structure. In addition, the codes can not be readily extended to 
include this interaction. Thus, the method does not take full advantage of these 
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specialized and well-trusted programs. The extensions of the multifunctional capability 
will focus on the loosely coupled method. 

The procedure for a loosely coupled method is given by (1) advance the structural 
system under the fluid-induced load, (2) transfer the motion on the wet boundary ( e.g 
the fluid-structure interface) of the structure to the fluid system, (3) update the fluid 
dynamic mesh accordingly, (4) advance the fluid system and compute new pressure and 
fluid stress fields, and (5) convert the pressure and stresses into structural loads. The 
multifunctional approach is applicable to steps two and five in the procedure outlined. 
These steps are concerned with the transfer of data from a CFD grid to a CSM grid. Data 
transfer is complicated by the fact that there are basic differences between the nature of 
the solution methods. CFD analyses are concerned with the flow field surrounding the 
surface exposed to the flow. Thus, a CFD grid is very fine around the exterior of an 
airfoil, wherever changes in the flow field characteristics ( i.e boundary layer effects) are 
expected to be maximum. Conversely, CSM methods examine airloads on the surface 
and how these loads affect the internal structure. CSM grids lie on the surface within the 
airfoil and are oriented to the structural components. Thus, CFD and CSM grids differ in 
grid density and data transfer requires extrapolation and interpolation of discipline- 
specific field variables. 

C'-y 

Smith et al. evaluated computational algorithms to interface CFD and CSM 
grids. In this reference, several candidate algorithms for passing information from the 
fluid regime to the structural regime were evaluated and the disadvantages of each were 
discussed. In addition, a load and motion transfer method based on the conservation of 
momentum and energy has been developed by Farhat 54 . In this reference, a conservative 
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algorithm for computing the loads induced by a fluid on a structure is discussed. This 
algorithm was shown to be accurate, robust and reliable for transferring data from a CFD 
grid to a CSM grid not only when the discretization differed, but also when the grids did 
not share the same geometry as in beam or wing-box geometric models (see Figure 5.19). 
In the figure, the structural surface is denoted by Ts and the fluid surface is denoted by 
Tf. The beam model is representative of the use of a beam finite element model to 
idealize the structural component within the airflow. The wing-box model is 
representative of a plate and shell finite element model to idealize the component in the 
flow. The multifunctional methodology developed herein provides an alternate 
conservative algorithm for transferring data from the CFD grid to a CSM grid. In 
general, the methodology can be used to transfer data among many different disciplines. 
Further development of the methodology to a two-dimensional (surface) interface is 
required. This development follows the approach presented by Aminpour et al. for 
coupling three-dimensional finite element meshes. 



Figure 5.19. Beam and Wing-box Structural Models. 

The governing equations for multifunctional analysis of vector-field problems 
have been developed in Chapter III and are given in Eqs. (3.32) and (3.34). Discretized 
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equations are given for solid mechanics in Eq. (3.57) and for fluid mechanics in Eqs. 
(3.60) and (3.61). In these systems of equations, the third equation represents the 
subdomain discretization mapping from one subdomain to another subdomain. This 
equation is given by 

K r « = 0 or [K^ K l2 H| = 0 (5.1) 

where the variables with a subscript 1 represent a solid subdomain and the variables with 
a subscript 2 represent a fluid subdomain. At this point, consider that the loads, a 2 , on 
the CFD grid are known. Eq. (5.1) can be used to solve for the unknown structural loads, 
« 1 , provided that matrix Ki, is square and invertible (i.e., the number of pseudo nodes 

used to describe the generalized displacement along the interface is equal to the number 
of Lagrange multipliers). Therefore, 

w l = - KjjK l2 « 2 (5.2) 

and 



- K i! K I 2 


^«2 ~ Aa 2 . 


Moreover, it can be shown that Kj A = 0 56 . That is, the matrix A spans the null space of 
matrix Ki. 

The fourth and fifth partitioned equations of the system of equations, given in 
Eqs. (3.57), (3.60) and (3.61), maybe used to interpolate the structural deformations to 
the fluid grid. Recall that these equations are associated with the generalized 
displacements on the interface and thus, the generalized displacement vector may be 


partitioned as 
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U 1 = 

■r 

► and similarly u 2 = < 

U 2 

/ 


L-IJ 


U 2 


where the subscripts, i and o, represent generalized deformations on the interface and 
within subdomain 1 or 2 ( e.g not on the interface). As such, the fourth and fifth 
equations are given as 

Kpj u j + uj =0 and Kp o u 2 + uj =0 


or 


Kpti*' +K i T u i =0 . 


Premultiplying this equation by A 1 yields 


A T K p u ; + A T Kj"u : =0 


Since A T Kf =0. 


A T K p u'' =0 and A T K pi uj +A r K po u^ =0 


or 


K iU {+K 2 i4=0 (5.3) 

The variables, uj , are associated with the known structural deformations from the 

structures grid, and the variables, u' 2 , are associated with the unknown deformations to 

be imposed on the fluid grid. Given that the matrix K 2 is square and invertible, Eq. 

(5.3) can be solved to obtain the unknown deformations. Therefore, 

u ; 2 =-k 2 'k,u; 


(5.4) 
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The values, u 2 , can now be used in the CFD code to update the surface deformation and 

to calculate a new set of surface loads. With Eqs. (5.2) and (5.4), the multifunctional 
methodology described herein may be extended to the multiple-domain analyses of 
different disciplines. 

5.5. SUMMARY 

In this chapter, the multifunctional methodology has been described and 
demonstrated for vector-field problems in engineering science. The selected problems 
included problems of solid mechanics and fluid mechanics. The governing equation in 
each case is the equation of linear momentum. In addition, for fluid mechanics, 
continuity conditions are required. The analyses performed have demonstrated the 
effectiveness and accuracy of the solutions obtained for the respective problems. In all 
cases, the results obtained using the multifunctional methodology were in overall good 
agreement with the reported analytical or reference solution. 

Based on the findings for the vector-field problems, extensions of the 
multifunctional collaborative methodology to multiple-domain analyses of different 
disciplines have been briefly investigated. An exploratory examination of the extensions 
illustrates the applicability of the methodology to loosely coupled multiple-discipline 


applications. 



202 


CHAPTER VI 

CONCLUSIONS AND RECOMMENDATIONS 

6.1. GENERAL 

Multifunctional methodologies and analysis procedures have been formulated for 
interfacing diverse domain idealizations including multi-fidelity modeling methods and 
multiple-discipline analysis methods. The methods, based on the method of weighted 
residuals, ensure accurate compatibility of primary and secondary variables across the 
domain interfaces. Methods have been developed for scalar-field and vector-field 
problems. The methods have been rigorously developed for multiple-domain 
applications, and the robustness and accuracy has been illustrated. Multi-fidelity 
modeling approaches have been developed that include both homogeneous ( i.e the same 
discretization method in each domain) and heterogeneous (i.e., different discretization 
methods in each domain) discretization approaches. Results have been presented for the 
scalar- and vector-field multifunctional formulation using representative test problems. 
Associated computational issues are also discussed. In addition, the extension to 
multiple-domain analysis with different disciplines has been discussed. 

6.2, CONCLUSIONS 

The multi-fidelity modeling of domains has been developed for homogeneous and 
heterogeneous discretization approaches for both scalar- and vector-field problems. The 
finite element and finite difference methods and combinations thereof have been used in 


each of the discretization approaches. 
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Multi-fidelity modeling 

Several general conclusions regarding the multi-fidelity modeling approaches can 
be made. First, each of the multiple-domain approaches leads to a non-positive definite 
system of equations, which impacts the solution strategy. Second, modeling flexibility in 
the multiple-domain method is increased at the expense of additional degrees of freedom 
introduced to the system of equations. However, the modeling advantage gained 
outweighs the computational expense due to the additional degrees of freedom, and the 
impact of the increased number of degrees of freedom due to the interface constraints is 
reduced as the overall problem size is increased. Third, while the multifunctional method 
encompasses heterogeneous discretization approaches using the finite difference method, 
the limitations regarding its use in the presence of complex boundary conditions and 
configurations restrict the method’s general-purpose use. Fourth, in general, the 
homogeneous and heterogeneous multiple-domain approaches using the finite difference 
discretization in one or both domains yield systems of equations that are not symmetric. 
This lack of symmetry is due to the use of the Dirac delta function as the weight function 
in the formulation. This function is introduced in the constraint integral used to form the 
coupling matrix in the upper triangular part of the system matrix. The finite difference 
“shape function” is used in the corresponding constraint integral used to form the 
coupling matrix in the lower triangular part of the system matrix. In fact, in the finite 
difference method, there may be a lack of symmetry in each of the independent 
subdomain “stiffness” matrices due to the imposition of the boundary conditions. 
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Scalar-field problems 

Conclusions regarding the multiple-domain modeling approach for the scalar-field 
problem include the following statements. First, scalar-field problems introduce many of 
the computational issues associated with the multifunctional approach. Second, 
satisfaction of the boundary conditions for the scalar-field problem using finite difference 
discretization is more straightforward than for the vector-field problem. The five-point 
template used to approximate the derivatives does not introduce difficulties at the corners 
of the domain, as is the case with the nine-point template used in the vector-field 
problem. Third, fictitious nodes are avoided by evaluating the governing equations only 
at the interior grid points of the domain. The essential and natural boundary conditions 
are applied at the boundary nodes with higher-order forward and backward difference 
approximations used for the first derivatives present in the natural boundary condition 
equations. Fourth, the governing equation is evaluated at the nodes along the subdomain 
common boundary. Straightforward central difference approximations are used at the 
interface to represent the interface tractions, which in turn are used to eliminate the 
fictitious nodes at the common boundary. 

Vector-field problems 

Based on the studies of the multiple-domain modeling approach for the vector- 
field problem, the following conclusions are drawn. First, the use of the finite difference 
method for the vector-field problem (e.g., plane stress problem) was far more 
complicated than for the scalar-field problem. The traction and displacement boundary 
conditions and the necessity to introduce and eliminate fictitious nodes outside the 
domain boundary greatly complicate the development. Second, the nine-point template 
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required in the finite difference approximation of the governing equations of the 
continuum introduces the need for alternative higher-order forward and backward 
difference approximations of the cross-derivatives present in the equations. Third, 
because of the difficulties associated with the first and second conclusions, the 
homogeneous and heterogeneous modeling approach using the finite difference method 
in one or both subdomains is not as attractive for vector-field problems as for scalar-field 
problems. Fourth, the governing equation is evaluated at the nodes along the subdomain 
common boundary. Complex manipulation of the nine-point template is required using 
forward and backward difference approximations of the cross-derivatives in order to limit 
the introduction of the fictitious nodes to the node along the common boundary at which 
the governing equation is being evaluated. This requirement is automatically satisfied in 
the scalar-field problem by the five-point template. The interface tractions are used to 
eliminate the fictitious nodes at the common boundary. 

Limitations 

While a rigorous multifunctional formulations has been presented, there are 
limitations in the implementation. Note that the purpose of the implementation described 
herein was to demonstrate the capabilities of the multifunctional approach on a set of 
representative benchmark problems. With this in mind, the limitations of the current 
implementation are as follows: 

• The nodes or grid points at the ends of the common subdomain boundary of each of 
the subdomains must coincide. 

• In the finite difference method used, at least three nodes are required in each of the 
coordinate directions where traction boundary conditions are imposed. 
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• Extreme care must be taken to perform accurate input and output using data-exchange 
fdes (in this work, double-precision floating-point accuracy). 

• The development of the interface routines in MATLAB limits the size of problem 
that may be analyzed. 

• Cubic splines are used on the subdomain common boundary, which requires at least 
four unique nodes along this boundary. 

• The implementation is limited to one-dimensional straight or curved common 
subdomain boundaries. 

• The geometry is assumed to be conforming. That is, each of the subdomains describe 
the same geometry along the common boundary. 

In this work, the benchmark vector-field problems illustrated require only C° 
continuity (continuity of the primary variable). Thus, continuity of the primary variable 
is maintained along the subdomain common boundary through the interface constraint. 
For plate bending problems using classical plate theory, C 1 continuity is required. In this 
case, continuity of the primary variable and its derivative is maintained along the 
common subdomain boundary. Here, the derivatives are approximated in the same 
manner as the primary variable. That is, cubic spline functions are used to approximate 
the generalized variables along the common subdomain boundary. Results for a wider 
range of problems including a plate bending problem have been given in reference 25. 
Summary of Results 

Results were presented for the scalar- and vector-field developments using 
example patch test problems. In addition, results for torsion, heat conduction and 
potential flow problems have been presented to demonstrate further the effectiveness of 
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the scalar-field development. Results for plane stress and plane flow problems have been 
presented for the vector-field development. Results for all problems presented are in 
overall good agreement with the exact or reference configuration by which they were 
evaluated. 

The multifunctional methodology presented provides an effective mechanism by 
which domains with diverse idealizations can be interfaced. This capability promises to 
provide rapidly the high-fidelity data needed in the early design phase. Moreover, the 
capability is applicable to the problems in the general field of engineering science and 
mechanics. Hence, the methodology provides a collaborative capability that accounts for 
discipline interactions among many disciplines. 

6.3. RECOMMENDATIONS FOR FUTURE WORK 

Future studies related to the present work are recommended. The present work 
provides a starting point for the following additional studies: 

1 . Explore the use of a finite difference energy method, which alleviates many 
of the issues associated with the proper identification of boundary conditions 
and the use of irregular grids. 

2. Evaluate the performance of the methodology for the analysis of more 
complex structures and fluid flow problems. 

3. Extend and implement the multiple-discipline capability. 

4. Develop other analysis capabilities including thermal analysis, modal and 
buckling analysis, dynamic analysis, and nonlinear analysis. 

5. Develop other heterogeneous multiple-domain discretization approaches such 
as the use of the finite element and boundary element methods. 
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6. Develop strategy to exploit massively parallel processing (MPP) computer 
systems. 

7. Incorporate computationally intelligent strategies to identify where and when 
homogeneous or heterogeneous approaches should be used. 
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APPENDIX A 

OVERVIEW OF STEPS IN ANALYSIS AND SIMULATION 

Multifunctional collaborative methods should address four typical steps of 
analysis and design, namely, (1) representation or modeling of the geometry, (2) 
knowledge-based selection and development of appropriate mathematical models (i.e., 
idealization/discretization), (3) solution of the mathematical model (continuous and/or 
discrete), and (4) interrogation/assessment of the results. These steps provide the 
foundation for enhanced integrated design and analysis tools, and the steps are briefly 
outlined in this appendix. 

Geometry Modeling 

To represent the structural geometry (geometry modeling) a geometric model is 
created to represent the size and shape of a system component. In aerodynamic and 
structural analyses, a common three-dimensional parameterized description of the 
airframe is shared. Geometry modeling is the starting point of the product design and 
manufacture process and is the first step in using a computer-aided design/computer- 
aided manufacturing (CAD/CAM) system ' The accuracy of the geometric model and 
the way in which it is structured has far-reaching effects on other CAD functions such as 
finite element analysis, drafting, and numerical control (NC) part programming. 
CAD/CAM systems can be utilized to develop a design and monitor and control the 

CQ 

manufacturing process from start to finish. Numerous CAD software packages for 
defining the geometry of structural systems are commercially available. 

Computer-aided engineering (CAE) has facilitated the assimilation of the 
engineer/analyst earlier in the design stage as an engineer in-the-loop. Typically, this 
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cycle leads from the design engineer to the analyst and back to the designer. A critical 
aspect of this cycle is the time required to generate analysis models, perform the analysis 
and decide if changes are needed. However, new trends in modeling and simulation are 
redefining the roles of the designer and the analyst. Many companies are now turning 
designers into analysts. The underlying philosophy guiding this paradigm shift is the 
desire to give designers the tools needed to predict a design’s performance early in the 
process, rather than just to define its geometry. These tools also embody a knowledge 
base to guide the designer through various analysis steps. Moreover, this new paradigm 
allows the highly specialized analysts to impact the design by performing more complex 
analyses to determine the structural integrity, the potential failure mechanisms and the 
complex response characteristics (i.e., material or geometric nonlinearity), and 
multidisciplinary characteristics of the design. 

This role redefinition can succeed only if enough analyses are performed early in 
the design process to identify critical design parameters, evaluate their interactions, and 
determine the best overall design. To expedite this process developers of computer-aided 
design (CAD) and analysis software have integrated the CAD and analysis functions. 
Such software integration and database coupling frequently enables designers to perform 
analyses directly on geometry, thus reducing the time required to prepare analysis 
models. 

Idealization/Discretization 

To develop discretized mathematical models of aerospace systems, several 
approximate numerical analysis methods have evolved over the years. The most 
commonly used discretization methods are the finite difference method and the finite 
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element method. The finite difference method of a configuration gives a pointwise 
approximation to the governing equations. While finite difference techniques are widely 
used in fluid dynamics and can treat fairly complex problems, they become hard to use 
when irregular geometrical shapes or unusual boundary conditions are encountered. This 
adverse attribute is particularly significant in structural analysis. In contrast, the finite 
element method is widely used for the analysis of many engineering problems involving 
static, dynamic and thermal stresses of structures. Typical input for a finite element 
analysis program consists of the geometric idealization, the material properties, the 
loading, and boundary conditions. The area of greater difficulty in the finite element 
technique lies in the geometric idealization, that is, representing the geometry of the 
structure by a suitable finite element mesh. Element aspect ratio, taper, and skew are 
characteristics that adversely affect the performance of many finite elements in use today 
and thus are factors in determining the suitability of a mesh. As the complexity of 
structural configurations and material systems being modeled with the finite element 
method has increased, manual mesh generation has become extremely tedious, time- 
consuming, expensive and consequently, intractable. This limitation is alleviated through 
the development of automatic mesh generators, which are typically integrated within the 
finite element modeling software and often integrated within the CAD system. These 
mesh generators are powerful tools for discretizing complex structural configurations. 
Issues associated with idealization still arise such as whether to use solid finite elements 
or shell finite elements. However, if the CAD and analysis engines are not driven from 
the same geometry, the translation of geometry may introduce errors in analytical models. 
In addition, due to the geometric complexity of such configurations, even the most robust 
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automatic mesh generator can often require analyst interaction to establish a suitable 
mesh and to provide engineering insight into the proper finite element to be used in the 
analysis. For example, some automatic mesh generators place three-dimensional models 
where two-dimensional shells should be used, which may distort the results. 

Response Prediction 

To solve the discrete system of simultaneous equations resulting from the 
discretization process and subsequent finite element assembly operations, myriad solution 
strategies have been developed for obtaining efficiently the unknown nodal values of the 
field variable or the primary unknowns. Two families of methods for solving linear 
systems of algebraic equations can be distinguished: direct and iterative equation solvers. 
The former can be defined as leading to the solution of a linear system in one step, while 
the latter will require many iterative steps. If the equations are linear, a number of 
standard solution techniques maybe used which generally include either an iterative or 
direct solver. If the equations are nonlinear, their solution is more difficult to obtain. All 
approaches will necessarily be repeated solution of linearized equations. A common 
solution method used to solve nonlinear systems of equations is the Newton-Raphson 
incremental-iterative solution procedure, which is accurate and converges for highly 
nonlinear behavior. High-performance equation solvers are a key component of solution 
strategies for linear and nonlinear structural response calculations for static, dynamic and 
eigenvalue problems in finite element analysis. There has been a plethora of research in 
the area of equation solvers for large-scale aerospace structures with only representative 
works referenced herein. Matrices resulting from discretization of structural systems are 
generally real, symmetric, positive definite, banded, and sparse. The performance of 
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iterative and direct equation solvers has been compared to identify the most appropriate 
tool for the solution of equations arising from structures systems 59 . This work identified 
advantages and disadvantages of both types of solvers. The study concluded that the 
relative performance of solvers depends on the amount of computations as well as the 
rate at which operations can be carried out on a given computer. 

Direct sparse solvers were found to be most attractive for models composed of 
higher-order finite elements, where they benefit most from a greatly reduced operation 
count. Sparse direct techniques are efficient improvements over first-generation direct 
methods that require more operation counts and larger memory capacity 60 . The number 
of operations in a sparse method are significantly reduced through reordering and storage 
strategies that effectively compress the global stiffness matrix into a format that exhibits 
a greater degree of nonsparsity prior to factorization and thus substantially reduces the 
associated computational costs. Iterative methods require much less memory than direct 
solvers, but their effective use depends on a fast convergence rate, which has been found 
to be best for finite elements with low aspect ratios. Skyline and variable band linear 
equation solvers have been developed to exploit the matrix characteristics of structural 
systems and to exploit the full capabilities of parallel and vector supercomputers 61 . More 
recently, general-purpose equation solvers have been developed for complex, 
nonsymmetric, indefinite, and dense matrix characteristics, which are prevalent in 
disciplines such as electromagnetic and acoustic analysis . Over the years, equation 
solvers have been developed to take advantage of the rapidly increasing computational 
power afforded by vector and parallel high-performance computers. These ultra-rapid 
equation solvers coupled with the major advances in computational power now available 
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in desktop personal computers and workstations have made it feasible to perform high- 
fidelity analyses in the preliminary design stage. However, additional developments are 
required to perform real-time large-scale analyses within an interactive virtual reality 
analysis and design environment. More intensive reviews of equation solvers may be 
found in the open literature ( e.g ., references 63, 64, and 65). 

Assessment of Results 

The fourth and final step in the analysis and simulation process is the 
interrogation of the results. In years past, the engineer would spend an enormous amount 
of time plowing through pages of computer output while waiting for results from 
additional analyses. With the increased speed and efficiency of today’s equation solvers, 
the rapid interrogation of results becomes decidedly more significant. It is at this step of 
interpretation of results that the engineer must be integrally involved. Powerful pre- and 
post-processing tools coupled with state-of-the-art computational technology provide the 
engineer with a comprehensive tool set for creating and discretizing complex geometries, 
performing analyses and visualizing results. Some software provide novel capability to 
enhance the designer-computer interaction while interrogating results. Engineers can 
view the results of parametric studies in a series of windows to identify or compare 
important design parameters. In addition, analysis results from different design 
approaches may be viewed in different windows and assessed to determine the most 
feasible design. This and other such visualization capabilities facilitate the rapid 
interpretation of analysis results, thus improving productivity of higher-order analyses 
and providing an opportunity for the engineer/analyst to be an integral part of the design 
process from concept to manufacture. Recently, immersive virtual reality environments 
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for visualization and interpretation of geographically dispersed results have been 
proposed as part of the NASA Intelligent Synthesis Environment (ISE) Initiative that 
promises to revolutionize the design process ’ . Immersive environments are human- 

scale computer-generated projection systems that allow users to interact directly with 

68 

their data in three spatial dimensions. Emerging advanced engineering environments 
will provide visual, auditory, and haptic feedback to further aid the engineer in detailed 
assessment of results. 
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APPENDIX B 

CUBIC SPLINE INTERPOLATION MATRICES 


The interpolation matrices used in the deformation and geometry assumptions of 
the multifunctional approach are outlined in this appendix. Given a series of points X/ 

(i = 0 , 1 ,... ,n) which are generally not evenly spaced, and the corresponding function 
values the cubic spline function denoted g(x) may be written as 
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where Ax=X/+ ; - Xj and g, XJt; denotes differentiation twice with respect to x. This equation 

provides the interpolating cubics over each interval for i = 0 , 1 ,...,/? - 1 and maybe given 
in matrix form as 


g = Tlg, xx +T 2 f (B.2) 

For each of the k values of x at which the spline function is to be evaluated, x/< x* <x/+ 7 , 

A /V 

k = 1,2 . . . p, and p is the number of evaluation points. The Tj and T 2 matrices may be 
written in the form 
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Note that there are, at most, two nonzero coefficients in each row of the T j and T 2 
matrices given above. 

Applying additional smoothness conditions (l.e., equating the first and second 
derivatives of adjacent interpolating cubics at Xj) yields a set of simultaneous equations of 
the form 


Ax,- 


i — 1 


A Xj 


Si xx ( x i —\ ) 


+ 


2 ( x i+l - x i) 
A Xj 


Sixx ('h ) Dk ’XX (*/+l) 


= 6 


f( x i+\)~f( x i) f( x i)-f( x i-l) 


i =\,2,...,n-\. (B.3) 


(A Xjf (Ax ( XAx m ) 

If the x/ are evenly separated with spacing Ax, then the Eq. (B.3) becomes 
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9joc ip^i— 1 ) ip^i ) "*" bk?xx (*^/+l ) 

J /(-V,tl)-2/(-V,)+/(-V,_i) l . (B.4) 

L (Ax) 2 

Eqs. (B.3) and (B.4) may be written as 

Ag, xc = Pf (B.5) 

The coefficients of matrices A and P are dependent upon the end conditions, which are 
discussed in the following section. 


End Conditions 

Whether the equations are of the form of Eq. (B.3) or Eq. (B.4), there are n- 1 
equations in the n+\ unknowns g, xx (xq ), g, xx (xj ), g, xx (x ri ). The two necessary 
additional equations are obtained by specifying conditions on g, xx (xq) and g, xx (x n ). For a 
natural spline, g, xx (xq)= g, xx (x n ) =0. However, in this work, these second derivatives are 
calculated by differentiating (twice) a cubic function which passes through the first four 
pseudo-nodes along the interface path and another cubic function that passes through the 
last four pseudo-nodes along the interface path. Evaluating this cubic function, 

g(x) = <3q + ct[X + ct2 x + a 3 x > ar| d at the first four points gives 

&( x o) 1 x 0 Xq Xq a {) 

g(*l) 1 X! xj 2 x\ < a\ 

g( x 2 ) 1 X2 xf x| a 2 

g{ x 3)_ 1 x 3 x 3 x 3 _ «3 

Solving for the coefficients yields a = N _1 g or 


or Na = g (B.6) 
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»11 n \ 2 »13 «14l[^o) 

«21 «22 «23 «24 < g(*l) 

*31 «32 «33 «34 g(x 2 ) 

»41 «42 «43 «44_|U( X 3) 

From the cubic function, g, xx (x) = 2a 2 + ba^x where a 2 and <23 are determined 
from Eq. (B.7). Equation (B.7) is valid for evenly spaced as well as arbitrarily spaced 
points. Similar expressions are obtained for the cubic function passing through the last 
four points where coefficients of the inverted matrix similar to those in Eq. (B.7) are 
denoted n \\ for k,l =1,...,4. With these end conditions, the matrices of Eq. (B.5) are 
given for equally-spaced points as 
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where p ^ = In^ + 6 / 74 ^ and p # = + 6 7j 4 £ for k,l =1,...,4. For unevenly spaced 

points, the tridiagonal A and P matrices may readily be obtained from Eq. (B.3). 
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Expressing g(x} in Terms of Functional Values fjxj X 

In Eq. (B.2), the spline function g(x) is expressed in terms of the functional values 
fix,) as well as second derivatives of the spline function, g, xx (Xj). However, it is desirable 
to express g(x) in terms of the function values f{xi) only. This manipulation is done by 
solving for g, xx (Xj) in Eq. (B.5) yielding 

g, xt = A" 1 Pf (B.8) 

Substituting in Eq. (B.2) yields 

g(x) = T 1 A _1 Pf + T 2 f = (fj A _1 P + T 2 )f = Tf . (B.9) 

Derivatives of the spline function are obtained by differentiating Eq. (B.9) yielding 
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0 for Xfr < Xf or > x, + j 

— — for Xj < x/ c < x j+ \ and j = i + 1 , 

Axi 

— — for Xj < x/ c < x i+ \ and /' = /' + 2 
Ax,- 

Again, note that there are, at most, two nonzero coefficients in the (f| ), v and 

matrices. In this derivation, x has been used as the independent variable. However, in 
the context of the interface definition herein, s is the independent variable and is 
substituted for x in the derivation in Appendix C. For the displacement assumption, the 
matrices developed for equally spaced points were used. For the geometry assumption, 
matrices for unequally spaced points were used. 
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APPENDIX C 

DERIVATION OF INTERFACE GEOMETRY 

C.l. GENERAL 

In the initial development outlined in reference 25, the interface path, P, was 
defined by piecewise linear segments. For curved interfaces, this definition only 
approximates the true curved geometry. The error in this approximation is a function of 
the interface path curvature and the number and location of the subdomain nodes along 
the interface. In addition, the interface path was computed along each subdomain 
independently, thus producing two different interface geometry definitions. For a 
structure with mild curvature, the error in the interface path definition did not influence 
the accuracy of the solution obtained in the analysis . However, for problems with 
moderate to large curvature, this error may be large and adversely influence the accuracy 
of the interface element analysis. 

In the present work, the element interface geometry is determined in one of two 
ways: (1) by specifying the function that represents the exact geometry of the interface 
(i.e., the linear interface is the trivial case) or (2) by passing a spline of the desired order 
(typically a cubic spline) through the specified coordinate data points to determine the 
function representing the geometry. In either case, the specified or computed function is 
parameterized and its first derivative is used to determine the arc length along the 
interface geometry of the subdomains as well as the interface boundary. Thus, in contrast 
to the earlier work, the interface geometry definition is a more accurate representation of 
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an arbitrarily curved geometry. In addition, only one interface path geometry is defined, 
and all the finite element nodes along that interface lie on that geometry. 

For a curved geometry, the most general way of determining the interface path of 
the two approaches mentioned previously is by using the latter approach (i.e., passing a 
cubic spline through the specified coordinates). In this case, a smooth curve is fit to the 
set of spatial coordinates by computing three cubic spline functions (one for each 
coordinate direction) expressing the coordinates as functions of a chordal distance 
parameter. The derivatives of these functions are obtained by differentiating the 
interpolating function. These derivatives are used in the parametric definition for the 
length of the arc between two points to compute the arc length between each of the 
specified coordinates. The spatial coordinates of the finite element nodes along each 
subdomain boundary provide the input for the interface geometry definition. These nodal 
coordinates are used to construct the function representing the curved geometry and to 
determine the arc length of the path. The associated variable, s, is computed along the 
subdomain boundaries. The number of evenly-spaced pseudo-nodes is determined 
internally or from the used-specified value after which the path variable, s, is computed 
along the interface path. See Appendix B for a brief discussion of the cubic spline used 
as the basis for the geometry representation. 

C.2. GEOMETRY REPRESENTATION 

The arc length or interface path is derived in this appendix. The spatial 
coordinates of finite element node i are given by x;, y ( -, and z/. The curve may be 


represented parametrically by 
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x = x(r) 

y=y (r) 

z = z(r) 

where r t = V ( x /+i _ ~i f + 0 ; /+ 1 -yif + ( z i + 1 - % f . Smooth cubic splines are fit 

through each of these coordinate functions. These coordinate functions are then 
expressed as 


x(r) = Tx s 
y(r) = Ty^ 
z(r) = Tz s 


where T is a matrix of interpolation functions (see Eq. B.9 in Appendix B) and is 
evaluated at the points r ( - The vectors x s , y s , and z s contain the sorted nodal coordinates, 
xuyu and z/, along the interface (i.e., the concatenation of the nodes from each of the 
subdomains to which the interface element is attached). 

The length of the arc between each of the points along the interface may be 
calculated immediately as 
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where T, r is obtained by differentiation of the interpolation matrix T with respect to the 
independent variable, r, (see Eq. B.10 in Appendix B) and is evaluated at points, r,\ The 
variable, s, is called the parameter of the arc length or the path variable herein. This 
variable measures the distance along the curve given by the parametric equations above. 
Thus, the arc length, s(rj) is obtained by numerical integration using Gaussian quadrature 
with four quadrature points. The path variable, as previously defined, is associated with 
the coordinates of the finite element nodes along the interface. The path variable, s, for 
the pseudo-nodes is computed by dividing the total arc length into equal segments. This 
total arc length is determined by summing the arc length between each set of two points, 
77 . i and /•;_ over the total interface path to obtain the total arc length. In addition to the 
path variable, s, at the j pseudo-nodes, the coordinate location of these pseudo-nodes is 
also desired. 

Moreover, in general, a computational coordinate frame is established along the 
interface; thus, the tangent to the interface path is desired. These calculations are 
addressed in the following discussion. 

Upon obtaining the path variable at the finite element nodes along the interface, 
the coordinate functions may now be expressed as 


x = x(^) = Tx s 

y=y(?) = Ty s 

z = z(s ) = Tz s , 

Here, the interpolation matrix T is evaluated at the path coordinates, s, of the pseudo- 
nodes yielding desired x, y, and z coordinates. The unit tangent vector to the interface 
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path is obtained by differentiating the coordinate functions with respect to the path 
variable, s, and is given by 

X ’S ( S ) ~ X, 

y*s (^) — 

z ’s ( s ) — Trv 2 s 

where T, s , is evaluated at the path coordinate, s, of the pseudo-nodes and the finite 
element nodes. The tangent vector is then given by 

j r< ( jV ) ? ! 3\v ( ,v ) J ^ z is{ s ) £ 

V r f ' r 

where > — (5)) + (}’•>$ (^)) + (^' )) . 
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